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ABSTRACT 


Several techniques for obtaining linear constant-coefficient airplane 
models from specified handling-quality time histories are discussed. One 
technique, the pseudodata method, solves the basic problem, yields specified 
eigenvalues, and accommodates state-variable transfer-function zero suppres- 
sion. The algebraic equations to be solved are bilinear, at worst. The dis- 
advantages are reduced generality and no assurance that the resulting model 
will be airplanelike in detail. 

The method is fully illustrated for a fourth-order stability-axis small- 
notion model with three lateral handling-quality time histories specified. 

The Fortran program which obtains and verifies the model is included and fully 


documented. 
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SYMBOLS 

a vector containing all of the elements of A 
A plant matrix with elements a^j 

b input distribution vector with elements b^ 

C coefficients matrix with elements c-j 

c final- value vector with elements c- 

«• I 

C* normalized longitudinal handling-quality criterion 

C*(t) longitudinal handling-quality criterion 

D( ) denominator polynomial 

D* normalized lateral handling-quality criterion 

D*(t) lateral handling-quality criterion 

d differential operator d 

G output measurement matrix, matrix transfer function 

6 augmented output measurement vector matrix 

h input/output coupling vector 

h augmented input/output coupling vector 

i pilot station to vehicle C.G. distance 

N{ ) numerator polynomials vector 

p normalized roll -rate handling-quality time history 

n 

p(t) roll rate 

• 

^co cross-over dynamic pressure in D definition 
r(t) yaw rate 

s Laplace transformation variable 

T similarity transformation matrix 

u(t) forcing function 
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V nominal airspeed 

state vector with elements x^(t) 
y output vector with elements y^(t) 

y specified response vector with elements y^(t) 

y augmented output vector 

3^ normalized sideslip handling-quality time history 

3(t) sideslip 

6, aileron deflection 

a 

6 elevator deflection 

e 

A eigenmatrix 

A specified eigenmatrix 

\ eigenvalue vector with elements 

<j)(t) roll angle 
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INTRODUCTION 

Aircraft control systems have been designed in the past to meet stability, 
apparent damping, and sensitivity criteria. Subsequently, pilot comments and 
ratings, using the Cooper-Harper rating scale, graded the performance of the 
pi lot-airplane-flight-control -system combination. The C-H scale emphasizes 

TAKEOFF 

CRUISE 

GROSS MANEUVERING AND AEROBATICS 

FORMATION FLIGHT 

TRACKING 

GROUND-CONTROLLED APPROACH 
LANDING APPROACH [1]* 

Experience has shown that the design criteria are often a subset of the han- 
dling-quality criteria. In order to design flight-control systems which 
achieve low (desirable) C-H ratings, some explicit consideration of handling 
qualities must be incorporated into the design process. The work begun by the 
Boeing Airplane Company and extended by the McDonnell -Douglas Company has re- 
sulted in quantitative specifications which are intended to insure “good" 
aircraft handling qualities. These specifications take the form of envelopes 
within which selected time-histories must fall, e.£. , normalized sideslip, 
normalized roll fate, and two blended quantities, C* and D*, containing 
the acceleration cues at the pilot station. 

An evolutionary process has provided the basic design of most aircraft 
control systems. It is often the case that the selection of numerical values 

★ — 

The numbers in square brackets refer to bibliography entries. 



for the many parameters in a control system is more difficult than the 
establishment of the control system structure and modes of operation. To 
overcome this "tuning" problem, a model -matching optimization technique can 
be implemented on a large digital computer. Such a program adjusts specified 
parameters of a simulated airplane-flight-control -system combination so as to 
minimize some measure of the dynamic differences between the closed-loop air- 
plane and a low-order model of a desirable prototype airplane. As a partic- 
ipant in the ASEE-NASA Summer Faculty Fellowship Program in the summers of 
1974 and 1975, the principal investigator, working in the Vehicle Dynamics 
and Control Division of the NASA Flight Research Center, Edwards, California, 
developed a method of translating handling-quality time-history specifica- 
tions into prototype aircraft models. These models are suitable for subse- 
quent use with an FRC model -matching program: CONOPT. 

The connection between specified time histories which fall within the 
established envelopes and numerical values for adjustable parameters onboard 
the aircraft is, in principle , established. The method is explained in the 
following sections and the lateral -motion case demonstrated. 

PROBLEM STATEMENT 

The NASA Flight Research Center computer program CONOPT, which resulted 
from the addition of the MIT Model Performance Index Design Program OPT to the 
in-house control system analysis program CONTROL, is a model-matching package. 
The user specifies the model in transfer-function form. The Model PI tech- 
nique [2] requires that the model transfer function, if it has any zeros, 
should have the same number of excess poles over zeros as the closed-loop 
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airplane control system. .If it has no zeros then the number of poles of the 
model transfer function should be equal to or less than the number of excess 
poles over zeros of the closed-loop airplane. Another condition on the model 
transfer function is that it must have reasonable eigenvalues. Thus the model 
is a linear, constant-coefficient ordinary differential equation in operator 
notation, 

D(d)x(t) = N(d)u(t) . (1) 

with loose bounds on the ' .''Ots of, 

0(x,) = 0 . (2) 

and constraints on the coefficients of N(d). In order to obtain the entire 
matrix transfer function in one operation the model will be represented in 
state variable form, 

dx = Ax + bu . (3) 

The problem is to determine an A and b combination which satisfies the loose 
bounds and constraints above and which describes a low-order prototype air- 
plane model with "good" handling qualities. FRC program CONTROL can be used 
to obtain transfer functions from A and ^ for use by CONOPT. 

The handling-quality time histories, y^(t), are linear combinations of 
state variables, 

y(t) = Gx(t) . (4) 

If all the y..(t) fall within their handling-quality time-history envelopes the 
model is appropriate for use with CONOPT. 
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The rigid-body equations of motion for sytranetrical airplanes slightly 
disturbed from straight and level flight can be separated into two sets of 
four simultaneous first-order linear constant-coefficient ordinary differen- 
tial equations [3]. One set describes the longitudinal motion, or motion in 
the plane of symmetry of a normally configured airplane. The other set de- 
scribes the lateral motion. In each case, 

dx = Ax + bu , (5) 

where A is 4 x 4 and b is 4 x 1 if there is only one input. 

Only one longitudinal handling-quality criterion has been established: 

★ ★ 

C^. Three lateral handling-quality criteria have been proposed: D^, normal- 

ized roll rate, and normalized sideslip. An envelope for the first derivative 
of each handling-quality time history has also been proposed [4]. These eight 
envelopes are plotted in Figures 1 through 4. These handling-quality time- 
history envelopes imply standard inputs. The C response results from a step 
change in elevator position, 6^, and the lateral handling-quality time re- 
sponses result from a step change in aileron position, 6^^. 

METHODS 

An obvious procedure for obtaining models with handling-quality time 
histories which fall within specified envelopes is to: 

1. Guess a model in terms of A and b 

2. Calculate x(t) 

3. Form y(t) 

4. Compare y(t) and y(t) to the appropriate envelopes 
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5. If no envelopes are violated, stop 

6. Otherwise alter A and/or b according to some 
strategy and return to Step 2. 

This "brute force" technique would consume large amounts of computer time and 
may not succeed. It would be necessary to define a cost I'unctional which 
measures all excursions outside of the envelopes and which is suitable for use 
in a numerical optimization scheme. Next a connection between the elements 
a. . and b- and the cost functional would have to be established. Finally, 

I j I 

there is no assurance that the attainable minima of this functional would be 
zero. The addition of hard and soft constraints arising from the loose bounds 
on the eigenvalues and specification of the pole-zero excess of the transfer 
functions complicates the numerical optimization problem. In the fourth-order 
lateral -motion case there would be: twelve variables to be manipulated, four 
soft constraint equations to be approximately satisfied, n hard constraint 
equations where n is the sum of the specified pole-zero excesses, and six en- 
velopes to be matched. 

One simplification which circumvents the need for a cost functional and 
optimization strategy is the use of direct search. If any solution which sat- 
isfies the constraints and the envelopes is as useful as any other, systematic 
or random direct search is easier to implement and no more wasteful of computer 
time than optimization. Random direct search has been used to obtain second- 
order longitudinal airplane transfer functions. In this simpler situation 
the a.. : i,j = l,2 are randomly selected. If the resulting eigenvalues fall 

* J 

within the region shown in Figure 5 [4], the b^ : i = 1 ,2 are randomly selected 
and c* and C* calculated. The ranges of the random choices are very 
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influential. Sample results are: 

1. 558 random sets of a., yielded 22 acceptable pairs of real 

^ J 

eigenvalues and six acceptable pairs of complex eigenvalues. 

2. For 4 having acceptable real eigenvalues 118 random pairs of. 
b^ were required to find one pair which yielded time-responses 
which match the C* and C* envelopes. 

3. For A having acceptable complex eigenvalues 320 random pairs 
of b^ were required to find one pair which satisfied the 
envelopes. 

4. If one envelope was matched (and the other violated) it was 
as likely to be the C envelope as the C envelope. 

The application of direct search to the fourth-order problem is currently 
being investigated. This is potentially useful for the development of models 
which satisfy the longitudinal criteria, C and C , but less so for the lateral 
case. In order to develop lateral motion models and obviate the use of iter- 
ative numerical techniques the problem must be reposed. 

If one begins with specific time responses which fall within the envelopes 
and seeks a model whose output responses closely resemble the specified 
handling-quality time histories the problem no longer requires Iteration. If 
the eigenvalues are specified or obtained from the input responses and then 
frozen and if all elements of A and b are unspecified rather than just those 
elements which are normally nonzero for most airplanes, the problem is well- 
posed. This procedure has three disadvantages: 

1. The partial degrees of freedom represented by loose bounds on 
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on the eigenvalues are lost. 

2. The envelopes are not directly used and may be violated by 
the results. 

3. The resulting A is not constrained to have the zero elements 
and kinematic elements usually found in such models of 
airplanes. 

Experience with lateral -motion examples has shown that A may well have ob- 
viously unrealistic elements and still produce realistic transfer functions. 
This is sufficient for present purposes. 

The reposed problem is the following: given y(t), h(a^j) and G(a^j), 

find A and b such that if 


and if 


then 


with 


dx = Ax + b6g^ , 
y = Gx + hfi^ , 

- -*• “ a 

/\ 

y ~ y in some best sense, 
A = A 


( 6 ) 

(7) 

( 8 ) 


Note that the elements of G and h are known functions of the unknown elements 
of A and that G is normally not square, there being fewer elements in y than 
in X. 

The first step is to obtain an analytical representation of the input 
time histories which one specifies in the form of a table of discrete values 
uniformly spaced in time. If the model is to be a fourth-order constant- 
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coefficient linear ordinary differential equation, its solution must have 
the form: 


y^(t) = c.^e^'* + ^ S4®' ''1 


( 9 ) 


+ c 
®i * 

These coefficients and eigenvalues, c^^ and Xj, are obtained by least-squared- 
error fitting the above form to the specified discrete values of y^. This is 
a two-step process. First the X. are calculated. The calculation is based 
upon first and second differences between adjacent y. values and is strongly 
influenced by their precision. If the values are represented by three signif- 


icant figures the resulting eigenvalues are drastically incorrect even when 
the data are contrived and there should be an exact solution. As the number 
of significant figures is increased the calculated eigenvalues more closely 
resemble the correct values. However, if data obtained from sketched time 
responses or imprecise tables are to be accommodated, it is impractical to 
calculate eigenvalues from such data. In such cases, the eigenvalues which 
the model is to have must be specified. In either case, the c.^ are straight- 
forwardly calculated and the resulting fits usually quite good. Thus one can 
represent the input time histories, y(t)» as 


y 




( 10 ) 


where the elements of c are the c. . 

An example set of discrete input data points and their fitted represent- 
ations is shown in Figures 6 through 8. Reasonable responses were sketched 
onto graphs of the lateral handling-quality envelopes and then converted to 
tabular form by estimating the values at half-second intervals. The eigen- 
values used in the curve-fitting process were specified. The values used were 
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\i = -2.4 
X2 = -.003 
X3 = -.25 - j2.0 
Xi, = -.25 + j2.0 

Once the specified discrete values of the y. are approximately represented 

y = Ce-* + c 


This equation can be Laplace transformed into 


Y(s) 


N(s) 
dTsT ' 


( 11 ) 


Assuming a unit-step input, the plant and output equations can also be Laplace 
transformed into 


sX(s) = AX(s) + b/s 
Y(s) = 6X(s) + h/s 
which are manipulated to obtain 

X(s) = (si - A)"‘b/s = ^(s)b/s 

Y(s) = G$(s)b/s + h/s . (12) 

If equation 12 is rearranged to match the form of equation 11, the right- 
hand sides of equations 11 and 12 can be equated, coefficient by coefficient, 
to produce n(nH-l) simultaneous algebraic equations where n is the order of 
the plant equation, equation 6, and m is the number of elements in the output 
vector, equation 7. For low-order models this is a satisfactory technique. 
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The algebraic equations are of degree <n. In the fourth-order lateral- 
motion example there are sixteen nonlinear equations containing, at worst, 
quadruple-cross -product terms. A limited amount of experience with these 
equations has shown them to be numerically sensitive. An effort is being 
made to develop a computer program which will solve the fourth-order example 
equations with various constraints on the elements and eigenvalues of A. 

While this may prove to be the most satisfactory approach to the problem, the 
numerical uncertainties are sufficient to motivate the development of an al- 
ternative technique which is computationally simpler. 

A less numerically-sensitive approach can be developed by combining equa 
tions 7 and 10. 

y = Gx + h6^ 

X = G'V - G"*h6, 

“ — •" “a 

X = G’^Ce^^ + 6"^c - 
dx = G’^CAe^^ - 6“*hd6, 

If these expressions are substituted into eauation 6 one obtains 

G“^CAe-* - G"‘hd6, = AG"^Ce-* + A6"^c - AG"^h6 + b6. 
or 

(G"^CA - AG"^C)e-^ = AG'^c + G"^hd6g + (b - AG'^h)^^ • 


For this expression to be true 


G"‘CA = AG’^C 


and 


(13a) 
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(AG’^h - b)6^ = AG‘‘c 

If T = G"^C, then 

A = TAT“' 
and 

b = AG‘^h - c) . 

It is assumed that 6g(t) s 1 and G’^dd^ is neglected at this point but must 
be introduced as an initial condition vector when computing model time re- 
sponses. A model obtained from equations 14 and 15 will have the specified 
eigenvalues and match the analytical representations of the input data exactly. 

Unfortunately the above results assume the existence of 6" whereas, in 
general, G is not square. One can substitute the right pseudoinverse and ob- 
tain similar results except that T will be singular and the model eigenvalues 
are no longer equal to the specified eigenvalues since A and A are no longer 
similar. If G is square and nonsingular, equations 14 and 15 hold. There are 
two cases of practical interest in which 6 is nonsingular. First, one can 
specify n = m, that is, let the number of time histories establish the order 
of the model. Secondly, it may be possible to adjoin created pseudo-time- 
histories to y until m = n. One must have the means to create these pseudo- 
time-histories. They must be linear combinations of the model state varia- 
bles and the resulting distribution matrix must be nonsingular. 

In the lateral -direction, small-motion example, the state variables are 
roll rate, yaw rate, sideslip angle and roll angle. The handling-quality 

'ff 

indicators are normalized roll rate, normalized sideslip angle and normalized D. 


{13b) 


(14) 
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The normalized roll-rate input data can be integrated to produce normalized 
roll -angle pseudodata. For the example, this was accomplished by exactly 
fitting a tenth-order polynomial to the eleven discrete roll -rate data points 
shown in Figure 6. This polynomial was integrated and evaluated at the same 
values of time to generate the discrete pseudodata points shown in Figure 9, 

The continuous curve in Figure 9 results from fitting an equation of the form 
of equation 9 to the discrete pseudodata points just as is done for the normal 
data. This allows the expansion of ^ to ^ and G to G and assures the pon- 
singularity of G. This, in turn, allows calculation of A and b using equa- 
tions 14 and 15. The same straightforward opportunity to augment y until G 
is square and nonsingular does not exist in the longitudinal small -motion case, 
although a designer could impose specifications on the motion, in addition to 
the Boeing C* criterion, and thereby create a nonsingular 6. 

COMPUTATIONAL PROCEDURES 

D* is a weighted combination of aircraft sideslip which is considered 
the principal low-speed handling-quality parameter and lateral acceleration 
at the pilot station which is the principal motion cue parameter at high 
speeds [4]. In terms of a stability-axis system representation of small 
perturbations about straight and level flight for a normally configured air- 
craft with the pilot station on the longitudinal principal axis the expression 
* . 

for D is 

D*(t) = V(B + r) + Ar + C3q^,Q3 (^6) 

where c^ is a dimensional constant defined in Table 1. In terms of the 
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Note that this assumes that the plant matrix is 


A = 


^11 «12 
®21 


With no constraints on the elements a... Thus A is specifically not con- 

* J 

strained to be of the form 


A = 


S L, 4 0 

Np N, 0 

% -■» Ye 
10 0 0 


( 19 ) 


similarly 


b = 


and not 


0 


Thus the output distribution matrix G for the lateral motion case is of the 
form 


G = 


ss 


. 1/p 
0 

Dp/Dss 


0 

0 

D./D 


r' ss 


0 

^ / 3ss 
V^SS 


0 

0 

V'^ss 


( 20 ) 


where 


y^(t) = p(t)/P55 = p^(t) 


y2(t) = 3(t)/eg5 = 6^{t) 
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ygCt) = D*(t)/Dg5 - D*(t) (21) 

and 


h = 


The augmented output distribution matrix is 


0 

0 

a 


( 22 ) 


6 = 


ss 


1/p 

0 

0 


Dn/Dcc 
p ss 


0 

0 

0 

V^ss 


0 

1/3 

0 


ss 


V“ss 


0 

0 

V“ss 


(23) 


which is nonsingular if / 0. The augmented output vector, y(t), is 


y^(t) = p„(t) 

(24) 

y4(t) = D^(t) 

Since n^ elements of A and the n elements of b are to be established equation 
14 is equivalent to n^ scalar equations with n^ unknowns and equation 15 is 
equivalent to n scalar equations with n unknowns. This is not significantly 
different from the Laplace transform method in which one obtains n(m+l) alge- 
braic equations without having had to create psuedodata. However, the degree 
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of the simultaneous equations to be solved in the Laplace transform method 
is n whereas the equations arising from the pseudodata method contain cross- 
product terms at worst. This dramatically Increases the liklihood of obtain- 
ing solutions by iterative numerical means. For the fourth-order lateral 
motion example, there are sixteen equations to be solved for the a^^. A 
Newton-Euler procedure [5] was employed in which an initial estimate of the 
solution vector a is iteratively improved. If equation 14 is rewritten 

f(a) = g (25) 

and P(a) is the Jacobian matrix associated with the equation 25 then 

This simple procedure has proved to be so successful that A = I can be used as 

-5 

the starting point and every element of “ ®|^) is reduced to 10 in three 

iterations typically. The elements of b are obtained from equation 15 without 
iteration. 

Before A and b can be calculated, equations of the form of equation 10 
must be fitted to the input data. The resulting coefficient arrays C and c 
appear in equations 14 and 15. Note that T = G in the pseudodata method. 

A set of eigenvalues can be obtained from each discretized input time 
history. If 

= y^(zAt) - y^(zAt) 


and 



XjAt 
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then 


e, = C,. 


'z 




z . r 

Va + C^. 


■ 14^4 


y^(zAt) 


(27) 


Where 


z = 1 ,2,... (number of discrete values y^(zAt)). 

The c.. and c. can be eliminated by linearly combining the e^- One obtains a 
1 j 1 

linear set of simultaneous algebraic equations 


D'y = d‘ 


(28) 


where the elements of D' are first differences of the table of y^(zAt) values 
and the elements of d' are second differences. Unfortunately j if the y^(zAt) 
values are obtained from time-response sketches or are imprecise for any reason 
the y. calculated from equation 28 are uselsss. Frequently the Pj obtained in 
such cases have negative values from which no X. can be calculated. 

In lieu of eigenvalues obtained from the input data, it has been neces- 
sary , in all practical calculations, to use specified eigenvalues. In such 

cases the calculation of C is based on minimizing 

11 

z=x ^ 

where the lateral handling-quality time histories typically span five seconds 
and are discretized by taking values every half second for a total of eleven 
values per time history. The initial value of y(t) is forced to be zero by 

As a check on the entire computational process once an A and b have been 


(29) 
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calculated, 6 and h can be calculated, x{t)obtained by integration and ^(t) 
calculated. This y(t), obtained by integration, is plotted with the fitted 
y(t). They should be identical. The integration method is described in 
Take 'ashi, et. , page 103 [6]. The solution matrix is represented by a 
series expansion containing p terms. The number of terms is specified by a 
recipe attributed to Paynter 


^ (nq)Pe"^ 0.001 (30) 

where n is the size of t.he A matrix and q is the largest element of AAt. 

The example fitted and integrated handling-quality time histories are 
plotted in Figures 10, 11, and 12. 
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ZERO SUPPRESSION 

The Laplace-transformed analytical representation of the input data, 
equation 11, has one fewer zero than poles in each element of the niatrix 
transfer function G(s) where 


Y(s) = G(s)/s 


or 


= S-X. 




S-A, 


S-X, 


s-Xt 


This expression can be rewritten, omitting the i subscript 

(c,x,t 02X2 + 03X3 + 0^1X4)5^ - (C3V4* W4^‘=’'1^3 

+ c^ X| Xg + C-| X-| X^ + c^Xi X^ + ^ ^ ^ 4 ^2^4 ^ ''I ^2^1 

+ C2 XiX2)s^ - ^‘^l^l V4‘^‘^3^1^3^4^ Vl V4‘*‘ Vl^3^4‘^‘^2^2V4 

+ Cj X 2 X|jX^ + C2^2^3^4 "*■ C] ^2^3 ^2^1 ^ 2^3 ^ ^ 3^1 ^2^3 

+ C-J XgX^ + C2X^X2X^ + c^X"! X2X^)s + X^ X2X2X^ (31) 

To reduce the number of zeros in the transfer function, the following con- 
straint equations must be incorporated into the least-square minimization, 
equation 29. To suppress one zero 


^ 1^1 ^ ^ 2^2 ^ ^ 3^3 ^ *" 4^4 “ ^ • ( 32 ) 

To suppress two zeros 


c^X^ + €2^2 + C3X3 + c^X^ = 0 

and 


'34 


i ! ] 


*“3^3^ + + CiX^Xj + CjX^X^ + CiX.jX^ + c^X^X^ 

^ ^ 2 ^ 2^3 ^ ^ 3 ^ 2^3 ^ 2 ^ 2^4 ^/|X2X^ ^ ^ 1^1 ^2 ^ 2^1 ^2 ” ^ 

If three zeros are suppressed, only one parameter remains to be established 
by the minimization and the resulting representation of the input time history 
is inadequate. 

For low-order models, zero suppression has two detrimental effects on the 
models obtained. The first is the worsening of the agreement between the dis- 
cretized input time histories and their analytical representations. For ex- 
ample, the unconstrained least-squared-error method yields a satisfactory fit 
of the rcll-rate input, as shown in Figure 6. The continuous fitted curve has 
approximately the same relationship to the envelope as does the discretized 
input data. If the LSE method is constrained by equation 32 to have one less 
zero, one obtains the fit shown in Figure 13. If one further constrains the 
LSE minimization by equation 33, another zero will be eliminated and, for the 
example data, the fit is degraded to that shown in Figure 14. The same se- 
quence is portrayed in Figures 7, 15, and 16 for the sideslip input. The 
Figure 7 results are unconstrained, one zero is suppressed in Figure 15 and 
two zeros are suppressed in Figure 16. Important features of the discrete 
inputs can be preserved by weighting the appropriate errors, e^, in the LSE 
method. This may make the constrained fits more useful. Even the uncon- 
strained cases can be altered. The risetime of the fitted roll rate can be 
reduced, for example. However, the overall fit cannot be improved by weighting. 

The second detrimental effect of zero suppression is the increased ten- 
dency of the A matrix to contain unrealistic values. The A and b arrays re- 
sulting from unconstrained fits of the input time histories are 
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-1.36 

-4.17 

-20.22 

-.02 

-1.12 

1.89 

12.43 

-.04 

0.22 

-.80 

-3.43 

3.02 

0.20 

1.04 

5.61 

0.00 


5.06 

1.59 

-.12 

1.09 


The Newton-Euler procedure for obtaining the a.j converged in three iterations 
to a maximum variation of a..,- ,ment from one iteration to the next of 
.00001 or less. This model was obtained from the unconstrained fits shown in 
Figures 6, 7, and 8. 

If one zero in the roll-rate transfer function and one zero in the side- 
slip transfer function are suppressed, the A and b arrays are altered to 


-3.49 

4.73 

97.14 

0.52 

0.96 

-3.42 

-22.48 

-.07 

-.19 

0.30 

4.00 

0.02 

-.32 

2.44 

15.42 

0.01 


0.00 

1.06 

0.00 

0.87 


The Newton-Euler procedure again required three iterations. This model was 
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obtained from the fits shown in Figures 13, 15, and 8. 

If two zeros in the roll -rate transfer function and two zeros in the 
sideslip transfer function are suppressed, the A and b arrays are 


-9.34 

-1.66 

403.88 

2.75 

3.51 

-2.88 

138.62 

-.86 

-.21 

-.04 

8.73 

0.06 

-2.10 

1.45 

96.10 

0.59 


and 

b = 

Three iterations were sufficient and the model was obtained from the data 
shown in Figures 14, 16, and 8. 

The number of terms in the expansion of the state transition matrix speci- 
fied by Paynter's recipe, equation 30, is a function of the largest absolute 
value contained in A and this number is frequently larger than the practical 
limit of about 30 when zeros are suppressed. 

The additional information required to calculate A and b from the discrete 
time histories shown in Figures 7-16 is given in Table 2. The specified 
eigenvalues were 


0.00 

1.45 

0.00 

0.87 


= -2.4 
= -.003 


^3 4 " -.25 + j2.0 
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AIRPLANE HOOEL HITH SPECIFIED TIME HISTORIES 

FLIGHT AWD VEHICLE PARAMETERS 
AIRSPEED 612*2 FT/SEC 

CG TO PILOT STATION LONGITUDINAL DISTANCE 22* 2<» FT 

DIMENSIONAL CONSTANT FOR DSTAR EQUATION -.3190 CUBIC»FEET/Le»SECONDS»SQUARED 

DYNAMIC PRESSURE 331,8 L0/FT-SQU ARED 

ROLLRATE NORMALIZATION FACTOR .500 

SIDESLIP NORMALIZATION FACTOR 10,000 

DSTAR NORMALIZATION FACTOR .010 


Jetstar Parameters and Flight Conditions 


Table 2 


The eigenvalues and vehicle and flight parameters are approximately those of 
a Lockheed Jetstar, a four-engined utility transport [7], at 20,000 feet al- 
titude and Mach = 0.6. It should be noted that the discrete input time his- 
tories were obtained fromp^^, and responses sketched by an FRC engineer 
[8] and do not refer to a particular airplane or flight condition. Finally, 
FRC program CONTROL was used to calculate transfer function coefficients from 
the above A and b arrays in combination with the appropriate distribution 
matrices, calculated from equation 22, to verify the suppression of the speci- 
fied zeros. 

VERIFICATION 

Flight test data obtained from the FRC Lockheed Jetstar [9] provides a 
test case for the above model -generation procedure. For small, lateral varia- 
tions about straight and level flight, the Jetstar can be represented by 


-2.353 

0.735 

-11.050 

0.000 

-.057 

0.358 

3.836 

0.000 

0.026 

-.999 

-.205 

0.053 

1.000 

0.054 

0.000 

0.000 


5.650 
0.031 
-.001 
0.000 

, *Body axis nondimensional stability derivative parameters used in place of 

stability axis values for verification purposes only. 

[ 
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1. 

0. 

0. 

0. 

0. 

0. 

1. 

0. 

0, 

0. 

0. 

1. 

14.75 

-7.044 

-146.0 

32.14 


h = 

The eigenvalues were 

= -.24045 
Xg = -.00310 
X3 ^ = -.25428 + j2. 06475 

CONTROL was used to calculate the handling-quality time histories for the 
period 0.0 to 5.0 seconds. The discrete input data extracted from the CONTROL 
output is given in Table 3. The flight and vehicle parameters are given in 
Table 2. The normalization factors were arbitrarily selected. 


0 . 

0 . 

0 . 

0 . 


The responses of the resulting model are shown in Figures 17, 18, and 20 


the 

discrete 

input values of 

Table 3 

arrays which 

define 

the model 

are 


-2.359 

0.771 

-10.939 

-.003 


-.053 

-.359 

3.849 

0.000 

— 

0.025 ■ 

-1.005 

-.201 

0,053 


0.906 

-.021 

.278 

0.003 
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5.818 


b = 


-.053 

-.215 


0.050 








1.0 

0. 

0. 

0. 


0. 

0. 

1. 

0. 

G = 

0. 

0. 

0. 

1 . 


14.35 

-11.13 

-143.5 

32.61 


h = 


0 . 

0 . 

0 . 

-.023 


These results are close to the original arrays. The differences are attrib- 
uted to the truncation of the CONTROL output values to two or three signifi- 
cant figures to approximate imprecise or sketched input data. 
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DISCUSSION 

The problem of obtaining an A and b from specified output time histories 
is one of nonlinear, noisefree identification. Five techniques for solving 
this problem have been suggested. The minimization of a cost functional which 
measures the differences between a trial solution and the handling-quality 
time-history envelopes would consume a large amount of computer time and there 
is no assurance that such a cost functional would be sufficiently well behaved 
to have a useful solution. Similarly, one could consume a large amount of 
computer time seeking solutions by random direct search. A graduate student 
is currently working on a variation of direct search in which the sensitivity 
of the time-history errors to changes in the A matrix elements is calculated. 
Then a set of incremental changes to the A elements can be obtained. 

The Laplace transformation method is also being pursued by a graduate 
student. In this form the problem is quite easily formulated but is ill- 
posed. It remains to be seen whether this method can advantageously incor- 
porate the loose bounds on the eigenvalues or not. It also has the disad- 
vantage of being a two-approximate-step process. One obtains y(t) from y(t) 
and then y(t) from y(t) except that the latter two will not coincide as they 
do in the pseudodata method. The Laplace transformation method and the sensi- 
tivity matrix method do have the advantage that they can be made to yield A 
matrices of the form of equation 19. 

The remaining approach, the pseudodata method has two advantages. It 
contains only one approximation step and it is niimerically efficient. The 
disadvantages are that it is somewhat less general and yields unconstrained A 
matrices. Experience has shown that the transfer functions resulting from 
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the unconstrained A matrices resemble those produced by A matrices of the 
fo'"m of equation 19. The pseudodata method results in a v/ell -posed set of 
bilinear algebraic equations which yield an A matrix having the specified 
eigenvalues. 

The pseudodata method is the only method which readily achieves zero 
suppression. In the other methods zero suppression contributes to their ill- 
posedness making useful solutions even more numerically difficult to obtain. 
The utility of any zero-suppressed solution is called into doubt by the 
detrimental effect suppression has on the LSE fit of the discretized input 
data {see Figures 13, 14s 15, and 16). Relieved of the need to suppress 
transfer-function zeros, one might prefer one of the other methods. 
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APPENDIX A 


Use of Program AANDB 
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INTENDED OPERATION 

1. Put in normalized discrete data 

2. Obtain normalized responses 

3. Put in normalization factors 

4. Obtain the normalizing distribution matrix 

5. Obtain the non-normal i zed airplane equations 
in state-space form. 

6. Use CONTROL to obtain transfer functions 

Steps 1 through 5 are performed by an example Fortran program AAND8 for 
the fourth-order small -lateral -motion case. The structure of ti^e program is 
shown in Figure A-1 . The subroutines perform the following tasks: 

MUGEN: Calculates eigenvalues from input time-history data. 

CASE!* (with entry points CASE2 and CASES): Calculates the 

required to fit c,e^^^ + c^e^^^ + c.,e^^*' + c.e^'*^ + Cr 

I <1 O H b 

to the input time-history data. 

LISTER: Prints information from COMMON on demand. 

FITTING: Fits a polynomial to roll-rate data and integrates it to 

produce roll-angle pseudodata. 

MODEL: Calculates the plant matrix. A, input distribution vector, 

b, the output distribution matrix, 6, and the input/output 
coupling vector, h^, 

RESPONS: Integrates the plant equations to produce comparison time 

histories. 

_ 

CASEl does not suppress any transfer function zeros, CASE2 suppresses one zero 
and CASE 3 suppresses two zeros. 
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GI’APHS: Plots envelopes and time histories and their first derivatives. 

INVR: Inverts matrices. 

The first fourteen input data cards are read by MAIN. The final two 
input data cards are read by subroutine MODEL. The input data cards contain 
the following information: 

1st data card : NN(20); 2012 

IF NN(1) = 1: LISTER will print: data, NN, MCASE, NPLOT, NHIST, 

input time histories 

NN(2) = 1: LISTER will print eigenvalues obtained from input 

time histories 

NN{3) = 1: LISTER will print eigenvalues specified by user 

NN(4) = 1: LISTER will print roll angle data obtained from 

roll rate 

NN(5) = 1: LISTER will print complex coefficient matrix C 

NN(6) = 1: LISTER will print summary of time-history curve fitting results 

NN(7} = 1 : not used 

NN(8) = 1: not used 

NN{9) = 1: LISTER will print envelopes for PN, BETAN, DSTAR 

& derivatives 

NN(IO) = 1: LISTER will print fitted curves for PN, BETAN, DSTAR 

& derivatives 

NN(ll) = 1: LISTER will print number of terms included in series 

for e^^^ 

NN(12) = 1: LISTER will print difference equation P matrix 

NN(13) = 1: LISTER will print difference equation g vector 

NN(14) = 1: LISTER will print responses obtained by integrating 

model equations 

NN(15) = 1: LISTER will print derivatives of responses obtained 

from P, q 
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NN(16) = 1: 
NN(17) = 1: 
NN(18) = 1: 
NN(19) = 1: 
NN(20) = 1; 


LISTER will print "PN PLOTTED" if NPLOT requests it 
LISTER will print "BETAN PLOTTED" if NPLOT requests it 
LISTER will print "DSTAR PLOTTED" if NPLOT requests it 
not used 
not used 


2nd - 12th data cards : TIME, PN, BETAN, DSTAR; 4F10.0 


0.0 

0.0 

0.0 

0.0 

0.5 

PN(.5) 

BETAN (.5) 

DSTAR(.5) 

1.0 

PN(1.) 

BETAN (1.) 

DSTAR (1.) 

1.5 

PN(1 .5) 

BETAN (1.5) 

DSTAR(1.5) 

2.0 

PN(2.) 

BETAN (2.) 

DSTAR(2.) 

2.5 

PN(2.5) 

BETAN{2.5) 

DSTAR(2.5) 

3.0 

PN(3.) 

BETAN (3.) 

DSTAR (3.) 

3.5 

PN(3.5) 

BETAN(3.5) 

DSTAR(3.5) 

4.0 

PN(4.) 

BETAN (4.) 

* DSTAR(4.) 

4.5 

PN(4.5) 

BETAN(4.5) 

; DSTAR(4.5) 

5.0 

PN(5.) 

"iETAN(5.) 

DSTAR(5.) 


13th data card ; MCASE, NPLOT, NHIST; 314 

If NHIST ^ 0 subroutine RESPONS will be called to calculate the 
model time histories. If NPLOT 0 subroutine graphs will be called 
to plot the model time histories according to: 
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NPLOT 

0 no plots 

1 PN and PNDOT plotted 

2 BETAN and BETANDOT plotted 

3 PN, PNDOT, BETAN, BETANDOT plotted 

4 PSTAR, DSTARDOT plotted 

5 BETAN, BETANDOT, DSTAR, DSTARDOT plotted 

6 PN, PNDOT, DSTAR, DSTARDOT plotted 

7 all plots 

8 no plots 


MCASE must be appropriate for: 


GO T0(l,2,3,4,5,6,7,8,9), MCASE 


1 P(s) has no zero(s) suppressed 

B(s) has no zero(s) suppressed 

2 P(s) has two zero(s) suppressed 

3{s) has one zero(s) suppressed 

3 P(s) has one zero(s) suppressed 

3(s) has two zero(s) suppressed 

4 P(s) has one zero(s) suppressed 

3(s) has one zero(s) suppressed 

5 p(s) has one zero(s) suppressed 

3(s) has no zero(s) suppressed 

6 p(s) has no zero(s) suppressed 

B(s) has one zero(s) suppressed 

7 P(s) has two zero(s) suppressed 

B(s) has no zero(s) suppressed 

8 p(s) has no zero{s) suppressed 

B(s) has two zero(s) suppressed 

9 P(s) has two zero(s) suppressed 

B(s) has two zero(s) suppressed 
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14th data card : specified eigenvalues; 8F10.0 

1-10 First eigenvalues which is real 
21-30 Second eigenvalue which is also real 
41-50 Real part of third eigenvalue 
51-60 Imaginary part of third eigenvalue 
61-70 Real part of fourth eigenvalue 
71-80 Imaginary part of fourth eigenvalue 

15th card ; velocity, length, c^, B 53 , D*^; 7F10.0 

VELOCITY: Nominal airspeed in ft/sec 

LENGTH: Centerline length from CG to pilot in ft 

* 

c^: Dimensional constant for D 

2 

‘^co* dynamic pressure in Ib/ft 

^ss* rate normalization factor 

^ss* Sideslip normalization factor 

k 

^ss’ O^TAR normalization factor 
16th data card : Newton-Euler parameters; 14, F20.0 

ITMAX, 14, Maximum number of iterations of Mewton-Euler algorithm 

EPSI, F20.0, when every unknown (the elements of the A matrix) changes 
by an amount smaller than EPM, the Newton-Euler Algorithm 
stops 

Note : 50,0.00001 seem to work well. 

INTERPRETATION OF THE PRINTED OUTPUT 
If NN = 20*1 all of the following output will be produced: 
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LISTER calls from the main program print: 

1. Date, list table NN, case number, plot request code, time response 
code and the input time history data. 

2. Real and imaginary parts of each MU(4) and El (4) obtained from the 
input time histories. The El are eigenvalues, and each MU is y^ = e°‘^^i 
where 0.5 is the pniform time interval between input time-history data points. 

3. Real and imaginary parts of the specified eigenvalues and associated 
MU values. These are the eigenvalues used in the least-square fitting process 
to obtain C^, not the quantities derived from the data. The derived values are 
presented only for comparison purposes. 

4. The PHJ or roll-rate psudodata generated from the PN or roll rate 
input time history. The time intervals are the same as for the original data. 

5. The TIME RESPONSE COEFFICIENTS matrix, C ' , is printed. This is a 
4x5 array of complex numbers: 

^ + c , where C' = [Cjc] , 

which analytically represents the input time histories. This y(t) is printed 
as: FITTED TIME RESPONSES by a call to LISTER from RESPONS. 

6. A sumriiary printout of MU, El, C and c 
LISTER calls from subroutine RESPONS print: 

1. The enyelopes, upper and lower boundaries, for P„ 3^* ^ p^, ^ 

(j 'k 

and D . Values are given for every 0.1 seconds. These values are stored 
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in DATA statement? in the main program. 

2. The curves which have been fitted to the input time-history data are 

printed for p^^, and D . The analytic expressions for the curves arc differ- 

d d d ^ 

entiated and tabulated also giving p^, and -gr D . 

3. The number of terms taken to calculate the truncated series used to 
represent e-^'^ is PAYMTERS RECIPE NUMBER. See CONTROL, Takahashi, al. , 
page 103 [6]. 

4. The difference equation parameters P and q in: 



5. The numerically integrated response time histories where = 

6. The first derivatives of the integrated responses are tabulated at 
0.1 second intervals. They are obtained from 

Lh. = + bu,, 

and 

LISTER calls from subrouLIre GRAPHS prints: 

1, If WPLOT is such that plots are requested, LISTER will print "PN 
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PLOTTED" If and p.^ plots have been generated, "BETAN PLOTTED" if 6^^ and 

,1 J ^ 

^ plots have been generated and "DSTAR PLOTTED" if D* and ^ D’ plots have 
been generated. The envelopes are automatically added to each plot. 

SAMPLE PLOTTED OUTPUT 

Plots are produced in three groups which can be requested individually or 
in any combination. These groups, the PN group, the BETAN group and the DSTAR 
group, consist of a title line, a lower graph and an upper graph. Each lower 
graph shows the upper envelope boundary, the lower envelope boundary, the 
analytical curve which has been fitted to the input data and the time history 
obtained by integrating the equations of motion. These four traces are repre- 
sented by continuous lines and are plotted versus time on the horizontal scale. 
The fitted and integrated lines should coincide. In addition, the lower graph 
contains eleven discrete symbols representing the input data. 

Each upper graph shows the upper envelope boundary, the lower envelope 
boundary, the first derivative of the fitted analytical curve and the first 
derivative of the integrated time response. The last two should coincide. All 
four curves are represented by continuous lines and are plotted versus time on 
the horizontal scale. A PN group sample is shown in Figure A-2, a BETAN group 
sample is shown in Figure A-3 and a DSTAR group is shown in Figure A-4. The 
example plots are unusual in that the original data was generated by a linear 
simulation program using a fourth-order model. Thus one should expect excel- 
lent agreement between the input data and the fitted and integrated results. 
Data obtained from other sources will not be matched as well. 
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EXPERIENCE WITH THE PROGRAM 

Several bounds on program flexibility exist. The first involves the 
generation of eigenvalues from input-time history data. Subroutine MUGEN does 
calculate a set of from each set of eleven points representing one input 
time history. These values are not correct or consistent even when the data 
is carefully contrived for best results. This is due to the dependence of the 
X^ calculation upon differences between adjacent (in time) data values. In 
anticipation of the program being normally used with data oboa'ined graphically, 
the values used in the example were held to three significant figures. This 
so degraded the X^. calculation that the values from the three input time his- 
•'ories varied from the known values used to generate the data by two orders of 
magnitude in the case of small eigenvalues and occasionally had the opposite 
sign. In general, the complex eigenvalues were more closely identified than 
the real eigenvalues. This was particularly true if the data precision was 
increased. In view of the anticipated character of the input data and the de- 
sirability of specifying the eigenvalues of the resulting model, no further 
development of MUGEN is planned. 

A second difficulty arises when zero suppression is specified. Subroutine 
CASE2 and CASES calculate coefficients, C, which yield transfer functions 
having the correct pole-zero excesses. The resulting expression will normally 
be a poor representation of the input time-history data but that is unavoidable. 
Unfortunately, the suppression of zeros also tends to ill-condition the nu- 
merical equations which must be solved to obtain A. In extreme cases this pre- 
vents convergence of the simple Newton-Euler-Raphson algorithm employed in 
MODEL. 
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Finally, there is a general difficulty that may interfere with the in- 
tegration of the model equations of motion in RESPONS. If the A matrix is 
ill-conditioned an unv/orkable number of terms in the series approximation for 
may be required. This is caused by the fixed At of 0.1 seconds which is 
built into RESPONS. A limit of 30 terms is imposed. Other limitations v/ill 
undoubtedly come to light as experience with the program accumulates. 


63 



APPENDIX B 

Program AANDB Listing 
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FROGRAt' AANOB ( INPUT, OUTPUTl MAI 

COMMON / NAMES/NN( 20 ) , PN { tl ) , RETAN (11 ) , DSTAR (11 ) , TPN ( 11 ) MAI 

1 , PNF (51) , RET ANF (5i),DSTARF(51),PNlF(51) , BET ANIF (51 > , DSTMA I 

2 ARIF (51) , PNC (51 ) »BETANC(5i ) , DSTA FC(51 ), PNIC ( 5 1 ) , BETA N IC ( 51» ,0 STAR IMA I 

3C (51), APN (11) , ABET AN (11) , ADST AR(ll) ,BPN{11 ) ,9B£TAN(11 ) .BOSTAR (11 ) ,MA I 
4PHIN(11) MAI 

COMMON /PARAM/MU, EIGEN, C , NCA SF. , MCA SE, NZ, NPLOT, MAI 

INHIST, PNFINAL,BNFINAL,0SFINAL,TIME(51) ,P(4,4»,G(4) MAI 

CCMMON /LIMITS/PNU(51> ,PNL (51)»B£TANU(51),BETANL(51) »DST ARUC51), MAI 
1DSTARL(51),PN1U(51),PN1L(51) ,BETAN1U(51) , BETA NIL (51) , C STARIU ( 51) * MAT 
20STAR1L(51) MAI 

COMMON /FINAL/ P SS, B ET A SS, OST A R SS, OP, DR, OBETA, OPHI, A ( 4,4) , Bt 4) ,ODEMAI 
ILIA MAI 

COMPLEX MU(4) , EIGEN (4) ,C (5, 4 ) MAI 

DATA PNU/,0, .333,. 66 7, 1.0, 1. 043, 1.064,1. 076,1. 08, 1,08,1. 078,1,0 76, MAI 

11. 074. 1.0 72, 1,07,1.0 69, 1.067, 1.065, 1.0 63 ,1,0 61 , 1,059, 1 .0 57 ,1,0 55, IMA I 
2. 053, 1.0 51, 1.05,1, 0^8,1,043,1. 044, 1.04 2, 1,04, 1.0 38, 1.0 36, 1.03 4, 1,0 MAI 

3 3 2,1.03, 1.029,1. 027, 1. 025 , 1. 0 23 , 1,0 21 ,1 , 019 , 1. 017 , 1, 0 1 5 , 1. 013, 1. 0 IHA I 

41,l.Ql,1.0 08 ,1.00 6,t.004,1.002,1.0/,PNL/Q.O,0.01,O.C27,0.065,0.10 7MAI 
5 ,0. 16,0, 227, 0. 3, 0.35,0 . 40 9, 0.45 , C.479, 0 . 50 7, , 5 3, 0 . 55 6, . 5 79 , 0 , 6 , 0 , 6MA I 

62.0. 53 9, 0,65 6, 0.677, 0. 692, 0. 706,0, 72 2,0. 73 3, 0 ,746, 0.762, 0,77 8, 0,79 MAI 

7. 0. 8 . 04 , 0. 816, 0, 829,0. 84, 0,®52, 0.861, 0.67 1,0. 86,0. 88 8,0, 896,0, 90 8, 0 MA I 
a ,915 ,0.9 2,0. 92 9,0. 93 8, 0 . 947,0. 95 6,0 .965, 0.974, 0.98 3, 0 . 99 2, 1. 0 / MAI 

DATA B ETA NU/ 21^2. 0, 2.1, 2. 2, 2, 3 , 2, 4,2. 5,2. 6 ,2. 7 ,2 , 8 ,2 . 9 ,3,0 ,3. 1 ,3,2 MA I 

1. 3. 3. 3. 4. 3. 5. 3. 6. 3. 7. 3. 8. 3. 9. 4. 0. 4. 1.4. 2. 4. 3. 4. 4. 4. 5. 4. 6. 4. 7,4. 8, 4MAI 

2.9.5. 0 /, BET A NL/ 2 1^-2 . 0,-2. 1,-2, 2, -2. 3, -2.4, -2. 5, -2. 6, -2. 7, -2, 8, -2, MAI 

39, -3. 0, -3. 1, -3, 2, -3. 3, -3. 4, -3. 5, - 3. 6, - 3. 7, -3. 8, -3. 9, -4. 0,-4. 1,-4. 2 MA I 

4,-4,3,— 4.4,— 4.5t-4.6,-4.7,-4»8?-4»9,-5.0/ MAI 

DATA D3T ARU/ 21^4. 0, 4. 2, 4. 4, 4. 5, 4. 8, 5. 0,5. 2,5. 4, 5. 6, 5. 8, 6. 0,6, 2, 6. 4MAI 

1.6.6. 6.8 .7.0 ,7. 2, 7. 4, 7.6, 7. 6,B. 0,8.2, 6.4, 8. 6, 8. 8, 9.0, 9.2, 9.4, 9. 5?9MAI 

2.8, 10.0/,OSTARL/2i’-4.0, -4. 2, -4.4, -4.6, -4. 8, -5. 0 ,- 5. 2,-5. 4.-5, MAI 

36 , ~5» 8, - 6. 0,~6o2,-6. 4, -5. 6, -6. 8,-7 .0,-7 .2 ,-7«4,— 7 .6,-7 .8,— 8. 0,-8. MAI 

4 2,— 8.4,— 8.6,— 8.8,- 9.0,— 9.2,— 9.4,— 9.6,— 9.8,'— 10.0/ MAI 

DATA PNl U/4, 0,4. 0,1. 68 6, 1.122, 0.992,0, 908, 0.343, 0,79 7,0, 724,0. 663 , MA I 

10.613, 0, 556, 0. 510,0. 4 6 4 , 0 . 4? 1 , 0 . 37 9 , 0 . 34 g , 0 . 3 1 8 , 0 . 2P 5 , 0 . 272, 0 .257, MAI 

20.234.0. 222, 0, 211,0. 199, 0, 192, 0. 188, 0.18 4, 0. 130, 0.176, 0. 172 ,0.158»MAI 
,3 0, 1 6 5, D. 1 69 1 0 .165, 0.165, 0.157’, 0.142, 0.119,0, 107, 0.C9 6,0. C8 C, 0 . G73 , M A I 

40. 069. 0. 061, 0 , 05 4, 0. 046. 0. ^46, 0. 046, 0. 04 6, 0.046/, PNl L/0. 306,0. 711, MAI 


ID 

20 

30 

40 

50 

60 

70 

e 0 

90 
ICQ 
110 
120 
130 
140 
150 
160 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
270 
2 80 
290 
300 
310 
320 
330 
340 
350 
360 
370 
380 



)Q«2559“o2?4?“#29f ?239**o3279'^ 



1C 


6-.223t-.205i-.133,-.148,-.129,- 

y-cllSi""* il6t~» 1149'‘» 1129-* ilO^** 
B-. 0 93,-,0 969 - .095 09 3t -oQ91,- 

9-. 079,-. 077, -.075,-. 074,-. 07^/ 
DATA 8ETAN1U/6 . 0 , 6.0 ,4 . 80, 4, 15, 
11. 78,1.&1,1« 49,1. 35,1.26,1.17,1 
2 — 6.0 ,- 5.0,— 4.60, — 4.15,— 3.72,— 3» 
3-1. 78, -1. 61,-1. i*9, -1.35, -1.26,- 
40/ 

DATA DSTARIU/12. 0,12. 0 ,9. 65,8. 3 
13. 90, 3. 55, 3,22, 2. 93, 2.70, 2.52, 2 
20ST ARIL/ -12, D , -12, 0 , - 9. 60 , - 8. 3 0 

3- 4.20, -3 .90,— 3. 56, -3. 22,-2.98,— 

4- 2. 0C,3Q»-2.0/ 

READ 300, NN 

READ 100, (TPMd) ,PN{I),BETAN(I) 
READ 200 ,MCASE,NPLOT,NHIST 
IF(NN<1> .EG. 1) CALL LISTERd) 
CALL MUGEM(PW:;HU, EIGEN) 

IF (NN? 2) . EG. 1) CALL LISTER(2) 
CALL MUGEN{9ETA N, MU , EIGEN) 

IF (NN ( 2) . EG. 1 ) CALL LISTEP(2) 
CALL HL‘GEN(DSTAR, MU, EIGEN) 

IF (NN(2) . EG. 1) CALL LI5TER(?) 
REan 4 30 , EIGEN 
MU ( 1)=CEXP(E IGEN(l) *^,5) 
MU(2)=CEXP(ETG£N(2)’^.5) 

MU(3 ) = CEXP(EIGrN<3)’ .5 ) 
MU(4)=C0MJG{HU(3) ) 

IF(MN(3) .EG.l) CALL LI3Tz??<3) 
CALL FITTING (PN, TPN, APN) 

DO 10 1=1,11 
PHIN (I > = 0 , 

00 10 J=l,ll 

PHIN( I) = PHIN( I) +0 PN ( J) ^ (I > 

IF(NN( a).EO.n CALL LISTER(4)^ 

GO TO (1, 2, 3, 4, 5, 6, 7, 8, 9) , MCA': - 
CALL CASEl ( MU , FM, EIGEN ,_C ( 1 , 1 ) ) 
CALL CA5F1 (HU , '’ET AN,E IGcN, 0(1, 


.32 3, - .319, - .30 4, -.289, - . 27^-, - , 2t>l , 
.125, -.118, -.125, -.123, -.12 1,-0 120, 

0 109, -.107, -.10 5,“.103,-.1Q2, - .100, 
.089,-. 088, -.086, -.084,-. 082, -.081, 


3.72,3.28,2.91, 2.61,2.3 7, 2.10,1. 95, 
.12,1.07,1. 0 5, 1.0 0,30’!./, 9ET ANIL/ 
28 ,-2. 91 ,-2 . 51 , -2 .37 , -2 . 10,-d9 5, 
1.17, -1.12, -1.0 7, ~i. 0 5,-1. 00,30 ♦-!. 

0, ?* 44,6. 56 ,5. 8? , 5.22 ,4.74,4.20 , 

.3 4, 2 , 24, 2. 14, 2 , 1 0,2. 0 0 ,3 0*2. 0 / , 
,-7.44,-6.5 6 ♦- 5. 6 2, -5. 22 ,-4.74, 

2, 7 0,-2. 52, -2. 3 4,-2.? i*, -2. 14, -2. 10, 


.□STAR (I) ,1 = 1, 11) 


** (12-J ) /FLO AT ( 12-J ) 


2 ) ) 


MA I 

MAI 

MAI 

MAI 

MAI 

HA I 

MAI 

Mi I 

MAI 

MA I 

MA I 

MAI 

MAI 

MA I 

MAI 

HA I 

MA I 

MA I 

MAI 

MAI 

MA I 

MAI 

HA I 

MAI 

MAI 

MAI 

MAI 

MA I 

MA I 

MAI 

MA I 

MA I 

MAI 

MAI 

MA I 

MA I 

MAI 

MAI 

MA I 

MAI 


3 90 
400 
410 
420 

4 30 
440 
450 
460 
470 
480 
49 0 
500 
510 
520 
530 
540 
550 
560 
570 
580 
590 
600 
610 
620 
530 
54 0 
650 
660 
670 

6 30 
590 
700 
710 
720 
73C 

7 40 
750 
760 
7 70 
780 



CALL CASE1CHU,PHIM» EI6EN» CI1»3H 
call CASEl CMU,0 STAR, eigen, C<1,4»I 
IF«NN(5».E(1. II CALL LISTER15I 
GO TO 11 

Z CALL CASE3CMU,PN, EIGEN, Cfl,lJI 
CALL CASE2(MU,BETAN, EIGEN, C(i,2»» 
CALL CASE1{MU,FHIN, EIGEN, CU,2II 
CALL CASEl(MU,OSTAR, EIGEN, C(1,4I> 
IF(NN<5) oEQ. II CALL LISTER(5I 
GO TO 11 

3 CALL CASE2(MU, PN, EIGEN , C(l,ll> 
CALL CASE3 (MU, BETAN, EIGEN, C(l,?n 
CALL CASE1CMU,PHIN,EI6EN, C(l,3n 
CALL CASE 1 (MU,CSTAR, EIGEN, C(l,4n 
IF(NN(5) ,£0.11 CALL LISTERC5) 

GO TO 11 

4 CALL CASE2(MU,PN, EIGEN, C(1,1M 
CALL CASE2(MU, BET AN, EIGEN, G(l,2»> 
CALL CASEKMU, PHIN, EIGEN , C(l,3n 
CALL CASEl (MU, DSTAR, EIGEN, C(l,4n 
IF (NN(5I ,EQo 1) CALL L1STER(5) 

GO TO 11 

5 CALL CASE2(MU, PN, EIGEN, C(l,l)l 
CALL CASEKMU, BETAN, EIGEN, C(i,2M 
CALL CASEl (MU, PHIN, EIGEN, C(l,3)) 
CALL CASEKMU, CSTAR, EIGEN, C(l,4)l 
IF(NN(5) .E0,1) CALL LISTEP(5I 

GO TO 11 

6 CALL CASEKMU, PN, EIGEN, C(l,in 
CALL CASE2(MU, BETAN, EIGEN, C(l,2n 
CALL CASEKMU, PHIN, EIGEN, C(l,3l» 
CALL CASEKMU, DSTAR, EIGEN, C(1,4I) 
IF(NN(5) • EQ, II CALL LISTER(5» 

GO TO 11 

7 CALL CASE3( MU, PN, EIGEN, C(l,ll> 
CALL CASEKMU, RET AM, EIGEN, C(1,2I) 
CALL CASEKMU , FHIN,E IGEN, C(l,3)l 
CALL CASEKMU, DSTAR, EIGEN, C«1,4II 
IF(NN(5» ,ECI, 11 CALL LISTER(5I 

GO TO 11 


MIX 

790 

NAI 

600 

N4I 

SiO 

Mil 

820 

MAI 

8 30 

HA I 

840 

MAI 

890 

MAI 

860 

MAI 

870 

MAI 

880 

MAI 

890 

MAI 

900 

MAI 

910 

MAI 

920 

MAI 

930 

MAI 

940 

MAI 

950 

MAI 

960 

MAI 

970 

MAI 

980 

MA I 

990 

MAI 

1000 

MAI 

1010 

MA I 

1020 

MAI 

1030 

MAI 

1040 

MAI 

1050 

MAI 

1060 

MAI 

1070 

MAI 

1U80 

MAI 

1090 

MAI 

1100 

MAI 

1110 

MA I 

1120 

MAI 

1130 

MAI 

1140 

MAI 

il50 

MA I 

1160 

MAI 

1170 

MAI 

118C 



8 CALL CASEH MU* PN, EIGEN, CJl.HJ 
CALL CASE3([MU*QETAN, EIGEN* CU,2>1) 

CALL CAS£i«MU*PHlW,EIGEN* C<1*3H 
CALL CASEU MU* OSTAR* EIGEN, C<l*^ll 

9 CALL CASE3tMU,PN,nGEN, 

CALL CASE3<MU,9ETAN, EIGEN, C«l,2n 
CALL CASEUMU,PHIN,EIGEN, CC1,3»? 

CALL CASElfMUeDSTAR, EIGEN, Cfl,<»M 
IFCNN(5I oEQolJ CALL LISTERC5) 

11 IF<NNf 6) oEQ.l) CALL LISTER«6) 

CALL MCOEL(C,EIGEN,PSS ,BETASS,tlSTARSS*OP,DR,OBETA,DPHI,A,B* 

IF(NHIST,NE, OJ pALL RESPONS 
IF{ NPLOT.NE. 0) CALL GRAPHS 
STOP 

100 FORMAH4F10, 0 » 

200 FORMAT (3 I4> 

300 FORMAT (2012) 

400 FORMAT(8F1Q,0) 

END 


MAI 
MAX 
MAI 
MAI 
M 1 
MAI 
MAI 
MAI 
MAI 
MAI 

OOELTAJMAI 
MAI 
MAI 
MAI 
HA I 
MAI 
MAI 
MAI 
MAI 


1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

13T0 



• 1 ^ 


@1 
^ 3 


so ^ 


SUBROUTINE MUGEN (0,HU, EIGEN) 

OIHENSION Dill) ,P?i>*4) •QI^»)»PI<'*»4JiiA«4» 
COMPLEX MUI4)«EIGENC4)«SflVE 
PC1»4)=D (2)-nU) 

P( 1 , 3 )=D ( 31 -D 12 ) 

PC1»2) =P(2,4) =D(4)-D«3) 

Pa,l)=PC2,3) = D(5)-D(4) 

P < 2 , 2 )=P ( 3 » 4 )=Q(l)=DC 6 )- 0 f 5 ) 

P{2*1)=P(3,3) =D(7)-Q (6) 

P( 3 , 2 )=P( 4 , 4 )=Q {Z) = 0 ( 8 )“ 0 ( 7 ) 

PC 3 , 1 )=P C 4 » 3 )=D( 9 )- 0 <S) 

P( 492 >=Q( 3 ) = DU 0 )-D{ 9 ) 

P( 4 » 1 )=D( 11 )-D( 10 ) 

D 12 = 0 ( 9 M- 3 .* ( 0 ( 11 ) - 0 ( 10 )) 

C 1 { 4 )=D 12 -D (11 ) 

CALL TNVR(P»Pl 94 , 0 » 4 ) 

□ 0 1 1 = 1,4 
A(I)= 0 . 

DO 1 J=l ,4 

1 ft(I)=ft(I)+PI(I, J) *Q( J) 

C 

c 

c QUADRATIC SYNTHETIC DIVISION 

C NUMERICAL CALCULLS 

C W, Eo MILNE PAGE 53 

0 
C 

X=t), 

Y= 0 « 

A 0 =lo 
Al=-A (1) 

A 2 =-A( 2 ) 

A3=-A(3) 

A 4 =-A( 4 ) 

10 B0=A0 

31 =A 1 -K*B 0 
92 =A 2 -X* 91 -Y *00 
e 3 =A 3 -X*B 2 -Y*Bl 
P 4 = A 4 -X*B 3 -V »02 


MUG 10 
HUG 20 
MUG 30 
HUG 40 
HUG 50 
HUG 60 
HUG 70 
MUG 80 
MUG 90 
HJG 100 
HUG 110 
MUG 120 
HJG 130 
HUG 140 
M J G 1 50 
MUG 160 
HUG 170 
HUG 180 
MUG 190 
HJG 200 
HUG 210 
HJG 220 
MUG 230 
HUG 240 
MUG 250 
MUG 260 
HUG 270 
HUG 280 
MUG 290 
HJG 300 
HUG 310 
HJG 320 
HUG 330 
HUG 340 
HJG 350 
MUG 360 
mug 370 
HUG 380 
MUG 390 



CO=BO 

C1=B1“X*C0 

C2=B2-X*C1»Y*C0 

C3= -X»C2-V*C1 

DD=C2*»2-C3*C1 

DX= (B3*C2“B4-*ClJ/00 

DV=(B^»*C2»B3*C3*/DD 

X=X*-DX 

Y=Y«-DY 

IF(ABStnX) «Ge,0, C001» GO TO 10 
IFIA9S (OY>.GE* 0. 0001 > GO TO 10 
U=B1 
V=B2 

ERR0R=A4-B2*Y 

IF(A8S(ERR0R1 .LE, 0. 0 001) GO TO 20 
PRINT 100, ERROR 
100 FORMAT(10X,6HERRCR=,F10.7) 

20 QUOT=X*»2-<+,»Y 
IFfQUOn 21,22,23 

23 ROOTl=-X/2,+SORT<QUOT)/2. 
R00T2=-X/2.-SQRT (QUOTl /2. 
NUtl»=CMPLX(ROCTl,0, 1 
MU(2)=CMPLX(ROOT2,Oa) 

GO TO 24 

22 MU(l»=CMPLX(-X/2*,0o » 

MU(2) = MU (U 
GO TO 24 

21 NU(l) = CMPLXt-X/2.,SQRT t-QUOT»/2. ) 
MU(2)=CMPLXi-X/2.,-SQRTf-QU0T>/2. I 

24 CU0T=U**2-4. *V 
IF(QUOT> 25,26,27 

27 R00T3=-U/2.*S0RT(QU0T) /2, 
RC0T4=-U/2o-SQRT (QUOTl /2o 
MU(3» = CMPLX(ROCT3,0, I 
MU(4)=CMPLX(ROCT4,0. I 
GO TO 28 

26 MU(3)=CMPLX(-U/2. ,0. ) 
h*U(4> = MlJ(3» 

GO TO 28 


BJIG 400 
KBG 410 
NJG 420 
((JG 430 
HUG 440 
HJG 450 
HUG 460 
HUG 470 
HUG 480 
HUG 490 
HJG 500 
MUG 510 
HJG 520 
HJG 530 
MUG 540 
HJG 550 
HUG 560 
MUG 570 
MUG 580 
HJG 5 90 
HJG 600 
HUG 610 
HJG 620 
HJG 630 
MUG 640 
MUG 650 
MUG 660 
HJG 670 
MUG 680 
HUG 690 
HJG 700 
HUG 710 
HJG 720 
HUG 730 
MUG 740 
MUG 750 
MUG 760 
HJG 770 
MUG 780 


;tn 



25 

MU«3»=CMPLXC»U/2«,SQRTI-QU0T»/2,> 

HJS 

790 


MU(4>=CHPLX!-U/2. «-SQRT J»QU0T»/2. > 

HJG 

800 

2B 

IF«A1HAG(HU(1H«EQ.0.) GO TO 30 

MUG 

810 


SAVE=HUflJ 

MJ6 

S20 


MU«1»=MU«3S 

HUG 

830 


MUC3»=SAVE 

HUG 

840 


SAVE-MU{2> 

HUG 

850 


HUC2)st1U(4» 

HUG 

860 


MUt4»=SAVE 

HJG 

870 

30 

CONTINUE 

HUG 

680 


DO 40 1=1,4 

HJG 

890 

40 

ErGEN(I) = CLOG«MU{I) ) ♦CMPLX (2. , 0 » ) 

HJG 

900 


RETURN 

HUG 

910 


END 

HJG 

920 



SUBROUTINE LISTER fWPRI NTJ LIS 

COMHON /NflMES/NNtSOJ ,PNlllUBETaNlliUOSTaR«liUTPM!lll LIS 

1 , PNF C 5 111 , e ET ANF < 5 U , OS T ARE 6 51 J * PNl F «5 i >, BET am F «5 1 J » OST LI S 

2aRlFf51> ,PNC tSlJ . 0 ETAMC « 51 » *0STARC«5i? ,PN 1 C «51l .8ETaNlC«51l*0ST«RlHS 
3C C51B, APN (11 1 » ABETANdl J » AOSTARf IIJ ,8PNI11 ) ♦BBETAN? IIJ .BDSTAR f 111 ,US 
4PHIWU11 LIS 

COMMON /PARAH/MUsEIGEN, C,NCaSE,MCASF, NZ, NPLOT, LIS 

1NHIST,PNFINALtBNPINAL9 0SFINAL# TIME (5111 ,P«4.A1 ♦ 0(A> LIS 

CCMMCN /LIMITS/PNUC51) ,PNL ( 5 U , BETANU ( 511 , BET ANL (51 1 , DSTARU (51 1* LIS 
lDSTaRL<51l9PNlU(51», PNIL (51 1 * BETANIUC 511 »BETAN SL(51) , DSTARIU ( 511 , LIS 
20STAR1L(51) LIS 

CCMMCN /FINAL/ PSS, BETA SS,DSTARSS *DP tOR.OB ETA, DPHI, A ( 4, 4) .Bt 41 » OOELI S 
ILTA LIS 

COMPLEX MU(41 ,EIGEN(4I ,C(5, 4) LIS 

PRINT 2000 LIS 

GO TO (10,20 *30,40*5 0, 60*70, 80, 90,100,110,120 ,130 ,14 0, 15 0,160, 170, LI S 
1180, 190,2001 ,NPRINT LIS 

10 CALL OATE(S) LIS 

PRINT 1000, S LIS 

PRINT 1010,NN,MCASE, NPLOT, NHTET,PN,TPN,BETAN, TPN,DSTAR, TPN LiS 
GO TO 210 LIS 

20 PRINT 1020, MU, EIGEN LIS 

GO TO 210 LIS 

30 PRINT 1030, MU, EIGEN LlS 

GO TO 210 LIS 

40 PRINT 10 40, PHIN, TPN LIS 

GO TO 210 LIS 

50 PRINT 1050, C LIS 

GO TO 210 LIS 

60 PRINT 1060, MU, eigen, C LIS 

GO TO 210 LIS 

70 CONTINUE LIS 

GO TO 210 LIS 

80 CONTINUE LiS 

GO TO 210 LIS 

90 PRINT 1090,(TIME(I>,PNU (I),PNL(I»,BETANU(n,BETANL(I» ,DSTARUII?, LIS 

IDSTARL(I) ,PN1U(I» ,PM1L(I» ,3ETAN1U(I» , BET AN 1L(I > ,OST ARiU ( I > ,BST ARl LLI S 
2CI),I=1,51» LIS 

GO TO 210 LIS 


ID 
20 
30 
40 
50 
60 
70 
80 
90 
100 
110 
120 
130 
140 
150 
1€0 
170 
180 
190 
200 
210 
220 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330 
340 
350 
360 
3 70 
380 
390 





00 


Q ^ 



100 PRINT 1100, CTINE«n *PNF8II»B€TAMF«IJ,0STARF«IJ ,PNlFeil ,BFTAIflF?II .LIS 


iiO 


120 


130 


140 


150 


10STARlF«n,I=l,51l 
GO TO 210 

1110, NZ 
210 

1120, t tP«I, JJ, J-1,4J, I=l,4» 

210 

11 30, a 
21 D 

1140, C TIME (II ,PNC!IJ ,eETaNC<I> ,OST ARC II H ,1=1,518 
210 

115 0,tTIHECI8 ,PN1C(I) ,BETAN1C(I) ,0STAR1C CI»,I=1,51I 
210 


160 


170 


180 


190 


200 

210 

1000 

1010 


PRINT 
GO TO 
PRINT 
GO TO 
PRINT 
GO TO 
PRINT 
GO TO 
PRINT 
GO TO 
PRINT 
GO TO 
PRINT 
GO TO 
PRINT 
GO TO 
CONTINUE 
GO TO 210 
CONTINUE 


1160 

210 

1170 

210 

1180 

210 


5HDfiTE=,AiO ,//J 


LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

US 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 

LIS 


,11F9,3,/,5X,7HTIME, ,11F9. 2 , //, 5X ,7H9ET AN, 

,11F9.2,//,5X,7H0STAR, , 11F9 , 3, / ,5X,7HT I ME , 


,11F9.3,LES 

,11F9,2LIS 

LIS 


RETURN 

FORMAT (I OX, 20 HLI ST PARAMETER TABLE, 201 3//5X12H CASE NUMBER 
1 PLOT C00E,I5,5X,19H TIME RESPONSE CODE, 15 , // , 20X ,1CHINPUT DATA, //LIS 
2,5X,7HPN , 

3/,5X,7HTIME, 

1020^FORMAT(5X,4*4X,4NREAL,10X,4HIHAG,11XI,//,2X,3HMU,,4(2E14.6,5X8,//,LIS 

103Q^FORMATC5Xa9HMU**VALUes°SPEC:FIEO,//,8X,4l 4HREAL ,1 OX ,4HIMAG, 10X8 LIS 

l,//,lX,3HMU,,4t2E14.6 8,//, IX, 6HEIGEN, ,4 IFIQ . 4,4X, FIO, 4,4X8 8 LIS 

1040 FORMAT (5 X,26HTHE GENERATED PHI DATA IS,, //,5X ,7HPHIN , ,11F9.3, 

i / 5 X 7 X ME 9 

1050 FORMAT85X,26HTIME RESPONSE COEFFICIENTS, // ,12X, 5 C 4X, ^REAL,6X,4HIMLIS 
1 AG, 5X8 ,//,5X,6HPN, ,5(2F10«4,3X8,//,5X,6HBETAN,,5C2F10 o4,3X8,//, L^S 

25X,6HPHIN, , 582F10. 4, 3X8, //,5X,6H0STAR,,5C2F10, 4,3X8 8 LIS 

1060 FCRHAT«5X,15HFITTING RESULTS, //,5X, 58 IIX ,4HREAL,6X ,4HIMAG8 


LIS 

LIS 


,//, LIS 


400 
410 
420 
430 
440 
450 
460 
470 
460 
490 
500 
510 
520 
530 
5 40 
550 
560 
570 
580 
590 
600 
610 
620 
630 
640 
650 
660 
670 
680 
690 
700 
710 
720 
730 
740 
750 
760 
770 
780 



15X,6MHU, *4f2F10,<»,5>O.//,5X,6HEIGEN,9<»I2F10,4,5XJ ,//,5X, LIS 

212HC0EFFICIENTS.//* 5X, 6HPN, * 51 2F10 • 4»5X8 , // 95X,6H9ETAN, , 5 1 2F10, LI S 

54*5X»,//,5X»6HPHIN* ,5 t2F10.4*5Xl «//,5X«6H0STAR, •5I2F10. 4«5X M LIS 
1070 F0RMflTi5X,28H«0DEL PflRAMETERS-SECCNO ROM.* 3F10.4, 5X, F10.4J LIS 

1080 F0RWflT(i5X,2rHM0DEL PARflfiETERS-THIRD ROHv^SFtO. «i*5X,F10,4» LIS 

1090 FORMAT IT5><:,18HRESP0NSE ENVELOPES ,/ i', 8X , AHTIME* 1 ?X »2HPN, 13X, 5H9ETAN, LI S 
113X,5HOSTAR,13X,5HPNOOT9lOX,SHBETANDOT,10X98HOSTAROOT,//951f3X«F9.LIS 
22,12F9.3*/>»/J LIS 

1100 F0RMAT<5X,21HFITTE0 TIME RESPONSES *//,8X,4HTIME 9 8X,2HPN, 5X, 5 HQETANH S 
l,5X,5HOSTAR,5X,5HPNDOT,2X,8HQeTANDOT,2X98HDSTAROOT,//»51i3X,F9,2,6LIS 
2F10,3*/>,/» LIS 

1110 FO RMAT(5X926HPAYNTERS RECIPE NUMBER LIS 

1120 FORMATf5X,29HOIFFERENCE EQUATION P MATRIX, »4F10. 3,/» 34X, AF10.3,/« LIS 
134X,4F10„3,/,34X,4F10,3» LIS 

1130 FORMAT (5X,29HDIFFERENCE EQUATION C VECTOR, , FlO , 3, /, 34X,F 10. 3 , /,34XLIS 
1,F10,3,/,34X,F10,3» LIS 

1140 FORMAT! 5X,19HINTEGRATE0 RESPONSE, // ,8X, 4HTIME.10X, 2HPN, 7X, 5H9-T ALIS 
1N,7X,5H0STAR,//,51 «3 X, F9. 2, 3F12.3 , / ) , /» LIS 

1150 F0RHAT!5X,41HFIRST DERIVATIVES OF INTEGRATED RES PONSES ,//,8X ,4 HTIMLI S 
1E,7X,5HPNDOT,4X,8H3ETANOOT,4X,8HDSTAROOT,//,51 !3X, F9 . 2, 3F12. 3,/» , ^LI S 


LIS 

IloO FORMATUQX,10HPN PLOTTED! LIS 
1170 FORMAT «10X,13HBETAN PLOTTED! LIS 
1180 FORMAT (10X,13HDSTAR PLOTTED! LIS 
2000 FORMAT!///! LIS 


750 
600 
810 
820 
8 30 
640 
650 
860 
870 
860 
890 
900 
910 
920 
930 
940 
950 
960 
970 
980 
990 
1000 
1010 
1020 
1030 
10 40 



03 

O 


8 | 

ft K) 
G ^ 


SUBROUTINE GRAPHS G?ft 

COHHON /NAHES/NN«20» ,PN!tl),BE1AN filUOSTARIllJ,TPMI li» GRA 

1 9 PNFC511 »BETANFI51> ,0STARFI51 I * PNlF iSl>,0ETANlF «51»«0STG?A 

2AR1F (51) »PNC C51!)*BFT ANC«515,0STARC<51I,PN1C«51 » ^SETANICJSJ. »• DSTARIGRA 
3C(51 KAPNtll )» ABETAMJll?#AOSTARgll» sBPNUlII •BSETANdlS »80STARmi *GRA 
ifPHINCllJ GRA 

COMMON /PARA M/MU, El GEN, C,NCASE,MCASE, NZ, NPLOT, 3?A 

1NHIST,PNFINAU,6NFINAL* 0SFINAL,T1MEI51> ,P<4,4> ,Q«4J GRA 

COMMON /LIMITS/PNUt5l» ,PNL«51> , 8ETANU «51 * , BETA NL C51 1 , 05TARU « 51 1 , GRA 
IDSTARL *51D,PNlU«51lvPNlL (51t ,BETANlU(!51II ,8ETANiLI515 ,0STAR1U(51» , GRA 
20STAR1U51J GRA 

COMMON /FINAL/ PSS, SET ASS,OSTARSS,DP,DR*OBETA, DPHI, A «4, A J ,B« , ODESR A 
iUTA GRA 

COMPLEX MU(4> ,eiGEN(4> ,C(5,4I GRA 

IF(NFL0T,EQ,1> GO TO 0 GRA 

IFCNFLCT.EQ. 3} GG TO IJ GRA 

IF(NPL0ToEQ,6J GO TO 10 GRA 

IF(NPL0T.NE,7> GO TO 20 GRA 

10 IF(NN(16)*EQ, II CALL LISTERtl6» GRA 

CALL QIKSET<5,0,0o0, 0.0, 3.0, 0,0,0. 5> GRA 

CALL QIKPLTCTIME,PNU,51,14H?TIME SECONDS?, 14HJPN 1/S ECONOSS, IIHSROSR A 
ILL PATEtl GRA 

CALL PLOTC-6, 0,1. 0,-3) GRA 

CALL QLIN£CTIH£,FNL,51,0) GRA 

CALL OLINE(TPN,PN,-ll, 74) GRA 

CALL QLINE(TIME,PNF,51.0) GRA 

CALL QLINE(TIM£, PNC, 51,0) GRA 

CALL PL0T(-1.5,3.0,-3) GRA 

CALL QIKSETfB. 0,0. 0,0. 0,3. 0,-2. 0,2.0) GRA 

CALL QIKPLT(TIME,PN1U,51,14H5TIME SECONDS? , BH$ OPN/DTS, 3HS $) GRA 

CA LL FL0T{-6. 0,1. 0,-3) GRA 

CALL QLINEtTIME,PNlL,51,0) GRA 

CALL QLINEfTIME,PNlF,51,0) GRA 

CALL 0LINE(T1ME,PN1C ,51,0) GRA 

CALL ENDPLT GRA 

20 CONTINUE GRA 

IF(NFLCT, EQ, 2 > GO TO 21 GRA 

IFtNPLOT.EQ. 3) GO TO 21 GRA 

IF(NFL0T.E0,5) GO TO 21 GRA 


10 
20 
30 
40 
50 
60 
70 
60 
90 
100 
110 
120 
130 
140 
15 0 
1€0 
17 0 
180 
190 
200 
210 
220 
230 
240 
250 
260 
2 70 
260 
2 90 
300 
310 
320 
330 
340 
350 
360 
370 
380 
390 



1F«NPLCT,NE,71) 6C TO 30 
21 .EQ* 1> CALL LISTER817J 

CALL QXKSET 8 5. 0 »O*0»O# 0 *3*0 *»5 .0«4*0) 


110H$SI0£SL1P$> 

CALL PLOT8“6« 0 9i*0««*3> 

(SLI NEJTIME9BETANL*51*0I 
QLINEaiHE»BETANF,51,0» 

QLINE(TPN ,BETAN»-11 «74> 
OLINE(TIMEteETANC,51,0» 
PLOTt-lo5,3«0, “3) 

QIKSETI5. OjO, 0,0. 0*3,0 »»6. 0,4. 01 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 


30 


PLOTC-6.0,1.0v- 3> 
0LlNEfTIM£*3ETANlL,51, 0» 
QLINE(TIHE,8ETAMF,51, OJ 
CALL QLINEtTIHE*3ETANlC,51,0» 

CALL ENDPLT 
CONTINUE 

TF(NPL0ToLT,4) GO TO 40 

IF( NN( 18) eSQ.l) CALL LISTER(IB) 

CALL QIKSET 85,0,0.0,0.0*4. 0*-10.0,5,0) 


113H$LATL. CRIT.Sl 
CALL PLOT (-6. 0*1.0, -3) 

QLINE8TIME,OSTARL*5l,C) 
0LINECTIME*0 ;TARF,51*0) 

QLlNEtTPN *r)STAR*-ll*74) 
OLIME8TIHE,OSTARC ,51*0 ) 
PLOTt-l.5,4.0,-3) 

QIKSETt5.0*0.0,0.0*3.0*-’i2.0*8,0) 


CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

1 ) 

CALL 

CALL 

CALL 

CALL 

CALL 


PL0TC-6o0,1.0*-3S 
QLINE«TIME,OSTAR1L,51*0) 
QLINE(TIME*0STAR1F,51* 0) 
QLINE( TIPE.OSTARIC ,51, 0) 
ENDPLT 


SR A 

%00 

GRA 

410 

6RA 

420 

S^A 

430 

G«A 

440 

GRA 

450 

G9A 

4E0 

6RA 

470 

GRA 

480 

GRA 

490 

GRA 

500 

GRA 

510 

iGRA 

520 

GRA 

530 

GRA 

540 

GRA 

550 

GRA 

560 

GRA 

570 

GRA 

580 

GRA 

590 

GRA 

600 

GRA 

610 

,GRA 

520 

GRA 

630 

GRA 

640 

GRA 

65 0 

GRA 

660 

GRA 

670 

GRA 

6 80 

GR A 

690 

GRA 

700 

$GRA 

710 

GRA 

720 

GRA 

730 

GRA 

740 

GRA 

750 

GRA 

760 

GRA 

770 


o cd o 
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O 


111 

;z 2: 

W 0^ 

^ 3 
3£ D 
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O 0^ UJ 


o 


82 
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SUBROUTINE CASei mu. 0* EIGEN, C9 

CAS 

10 


DIMENSION 0 811) .UCS .8} .UZNV 88,8 I.CC 881 ,PP8 8> 

CAS 

29 


DIMENSION UU86,6»,UUINV86,6),PPPPf6»,CCCf6l 

CIS 

30 


DIMENSION UUU84,4),UUUINV84,4i ,PPPPPP84) .CCCC84) 

CAS 

40 


COMPLEK MU841» ,D0 810 » ,Q 84,4 1 ,P(4 »,PC 819 ,4 f ,C 85 1 ,EIGEN 84) 

CAS 

50 


COMPLEX QQ83,3»,PPP 63J,RRC810,3I 

CIS 

GO 


COMPLEX QQat2,2J, PPPPP82> ,RRRCtiO,2l 

CAS 

70 


COMPLEX Al,fl2,A3,A4 

CAS 

ao 


DO 1 1=1,10 

CIS 

90 


DO 1 J=l,4 

CAS 

100 

1 

RC8I,J)= 8MU8J»<*I«1. 1 

CIS 

110 


DO 2 1=1,10 

CAS 

120 

2 

ooan=CMPLXf 8D8i+ii»ocin,o,» 

CAS 

130 


DO 3 1=1,4 

CAS 

140 


00 3 J=l,4 

CAS 

ISO 


D8I , J>=0, 

CAS 

160 


00 3 K=l,10 

CAS 

170 

3 

Q8I, J»=0(I, J»«^RC(K, J)*RC(K,II 

CAS 

ISO 


DO 4 1=1,4 

CIS 

19 0 


P(I»=0, 

CAS 

200 


00 4 K=1,1D 

CAS 

210 

4 

P(I )=P8I» +00 8K »*RC (K,I ) 

CAS 

220 


00 5 1=1,4 

CAS 

230 


DO 5 J=l,4 

CAS 

240 


Ua,J)=REAL(Qa,J)l 

CAS 

259 


U8 n-4 , J+4 J=REAL!Q{ I, JH 

CAS 

260 


U( 1+4 ,J»=AIMAG(08I,J» > 

CAS 

270 

5 

U8I, J+4 » = -AIMAG(QCI, JH 

CAS 

280 


CALL lNVRCU,UINV,e,0,8l 

CAS 

290 


DO 6 1=1,4 

CAS 

300 


PP tI) = REAL{Pf IH 

CAS 

310 

6 

PP(I+4l=AIMAG(PtIlll 

CIS 

320 


00 7 1=1 ,8 

CAS 

330 


CC(I) = 0. 

CAS 

340 


DO 7 J=l,8 

CAS 

350 

7 

CCCI»=CC(1>+UINV CI,J»*PP(J> 

CIS 

360 


no a 1=1,4 

CAS 

370 

8 

C( I>=CMPLXrCCC I> ,CC8 1+4 1 > 

CAS 

380 


C85»=-C(l)-C(2»-C(3»-C 841 

CIS 

390 



1 

oo 

-5:a 



RETURN 
entry CftSE? 
ft2=“EIGENfi2»/EI6ENUJ 
fl3=*EIG£N(I31l/EI6EN(tll 
A«»= = EIGEW848 /EIGEMf 
DO 10 1=1,10 

RRC«I,1» =1»*A2-NU (T21) *• 1° A2*HU« 1 > *•! 
RRC{I,2>=1.*A3-MU(3» »*I -A3*HU II J**I 
10 RRCII,31=1«<-A‘»-MU (4i** I-A4*MUI1»**I 
DO 20 I*lslO 

20 DDI IJ=CMPLXIDI I+lt ,0 .J 
00 30 1=1,3 
00 30 J=l,3 
QQII, J)=Q« 

DO 30 K=l,10 

30 00(1, J)=QOII» J>-RRCI K, J)*RRCIK,I» 

00 40 1=1,3 
PPP(I>=0, 

00 40 K=l,10 

40 PPP (I)=PPPCI 1 ♦C0CK)*RRC (K,I » 

DO 50 1=1,3 
00 50 J=l,3 
UU(T,J) = REALIQQ (I,JJ » 

UU(I +3,J + 3> = R£flL(QC(I, J) » 

UU( I-t-3 , J) =ATMAG ( OQ(I , J > > 

50 UU C I ,J+3 ) =- ftIMAG (00 1 1, J ) ) 

CALL INVR (UU,UUINV,6, 0, 6> 

DO 60 1=1,3 
PPPP (I) = REAL IPPPII* J 
50 pppp(l*3 J=AIMAG(PPPII) 1» 

00 70 1=1,6 
CCC (I)=0. 

DO 70 J=l,6 

70 CCC(II=CCC(I)+UUINV II,JI»PPPPCJ> 

DO an 1=1,3 

90 C(I«-lJ=CMPLX(CCCf II , CCC ( 1+3 H 

C (1) = A2*C (2)+A3*C(3I*A4*C (4) 

C (5 l=-C(l l-c (2 ) «*C (3 I -C (41 
RETURN 


CAS 400 
CIS 410 
CAS 420 
CAS 430 
CAS 440 
CAS 458 
CAS 460 
CA S 47Q 
CAS 480 
CIS 490 
CAS 500 
CIS 510 
CIS 520 
CIS 530 
CAS 540 
CAS 550 
CAS 560 
CIS 570 
CIS 500 
CIS 590 
CAS 600 
CAS 610 
CAS 620 
CIS 630 
CAS 640 
CAS 650 
CIS 660 
CIS 670 
CIS 680 
CAS 690 
CAS 700 
CAS 710 
CAS 720 
CfiS 730 
CAS 740 
CAS 750 
CIS 760 
CAS 770 
CIS 700 



ENTRY CASE3 

Al=EIGENl3»*JEIGEN«3l»EIGEN(l2>»/«EIGeN!ll*IEI6EN«2l»ei6EN«lHI 
A2=EIGENf4»*!EIGEN«4J-EIGEN«2? >/<EIGEM81l» «EIGEN!2»-ElGEtmil I 
A3= EIGEN 835* t€I6EN«3l •EIGEMfi 5 5 / 8EIGEN «2I* IEIGEN8 II “EIGEM82I 5 • 
fl4sEIGEN845*(EIGEN<<*5»EIGENCin/8EIGENf25*lEIGEN«15-EIGENC2lll 
DO 100 1=1,10 

RRRCa,15=l.*A3+ft 1-NUJ35**I»A3*HU<2S**I-A1 *HU«15**I 
IQO RRRCtl ,2l = lo<-A4*A2-MUf 45**1-“ A4*HU 825 ••1»A2*«U 815 
DO 200 1=1,10 

200 DO 8I5=CMPLX I081'»-15*0o 5 
DO 300 1 = 1,2 
DO 300 J=l,2 
QQQ8I,J5=0. 

DO 300 K=l,10 

300 QQQ8I, Jl =QQQ8I,J5“-RRR CK,J5*RRRC8K, 15 
DO 400 1=1,2 
PPPPP8I5=0, 

CO 400 K=l,10 

40 0 PPPPP8I5 = PPPPP (I 5*OD CK5*RRRC CK,I 5 

C 835 =8 000 82, 25 *PPPPP Cl 5 - QQQ 81 ,2 5 *FPFPP 82 5 5 / 8QQQ 81 *1 5 *QQQ 82,2 5- 
100081,25 *000 (2,15 5 

C 845= (QQQ81, 1 5 *PPPPP825 -0008 2, 1 5 *PPPPP 815 5 / 8000 8 1, 1 5 *Q0Q (2,2 5- 
1000( 1, 25 *0008 2,15 5 
C(15 = Al*C 835 +A2*CC45 
C 8 25 =A3*C 835 +A4*C (45 
C85 5 =-C8 15-C( 21-C 835 -C 84 5 
RETURN 
END 


CAS 

790 

CAS 

SOD 

CAS 

810 

CAS 

020 

CAS 

830 

CAS 

840 

CAS 

850 

CIS 

660 

CAS 

870 

CAS 

680 

CIS 

890 

CAS 

980 

CAS 

910 

CAS 

920 

CAS 

930 

CAS 

940 

CAS 

950 

Cl s 

960 

CAS 

970 

CAS 

980 

CAS 

990 

CAS 

1000 

CAS 

1010 

CAS 

1020 

CAS 

1030 

CAS 

1040 

GAS 

1050 

CAS 

10 eo 


00 
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SUBROUTINE INVR !A,8* JJJ « IT, MX J 
C lMX#MV,mJ,MS,MaTl,HaT2,HftT3, WftT<,,MflT5,MAT6l 
C PROGRAM AUTHORS R,E. FUNDERIC AMC R.6. EOMAROS, 

C COMPUTINF TECHNOLOGY CENTER, UNrC^ CARBIDE CORP,, NUCLEAR OIV*, 
C OAK RIDGE, TENN, 

C 

C CTC ORO PROGRAM NO, 9067.1 

DIMENSION A{HX,MX>,9<MX,MXI 

MAT1=MX 

MAT2=MX 

IF fJJJ.NE.lS GO TO 50 
0(1,U-1./A*l,11l 
RETURN 
50 CONTINUE 

00 21 1=1, JJJ 
DO 20 J=1,JJJ 
BCI, J)=OoO 

20 CCNTINIE 

Ba,i)=i.o 

21 CONTINUE 
KK=JJJ 
NV=JJJ 
0 = 1 , 

IF CJJJoLT. 0)0=0, 

KKM=KK-1 
00 9 1=1 ,KKM 
S=0,0 

00 1 J=I,KK 
r=ABS(A< J,I) ) 

IF CR.LT,S) GO TO 1 

S=R 

L=J 

1 CONTINUE 

IF (UoEQ,!) GO TO 5 
00 2 J=I,KK 
S=ACI,J) 

A(I, J) =A (L, J) 

4 <L, J) = S 

2 CONTINUE 


Itiv 10 

INV 20 
INV 30 
INV 4*0 
INV 50 
!<<V 60 

INV 70 
I1V 80 
INV 90 
IMV 100 
I«V 110 
INV 120 
INV 130 
I«V 140 
INV 15C 
INV 160 
INV 170 
INV 180 
INV 190 
INV 200 
INV 210 
INV 220 
INV 230 
INV 240 
INV 250 
INV 260 
INV 270 
INV 280 
INV 290 
INV 300 
INV 310 
INV 320 
INV 330 
INV 340 
INV 350 
INV 360 
INV 370 
INV 380 
INV 390 


Bp 


IF(NV.LE. OSGO TO 4 
DO 3 J=1,NV 
S“8«I* JJ 

B<L,JJ‘S 

3 CONTINUE 

4 D =-0 

5 IFCAfI,I».EQ,fl,»GO TO 9 
IPO=I+l 

DO 8 J=IPO,KK 
IF fAfJtIKEQ.O.) GO TO 8 
S=A(J9 I> 

A f J,It -0,0 
DO 6 K=IPO,K?( 

AtJ,KI=A( J,K»-A<I,K) *S 

6 CONTINUE 

IF (NV.LE.O) GO TO 8 
00 7 KslvNV 

8(J,K»=B (J,K)-3(I,K> *S 

7 CONTINUE 

8 CONTINUE 

9 CONTINUE 

DO 10 1=1, KK 
0=D*Aa,I) 

10 CONTINUE 
IF(NV,LEoO»GO TO 13 
KMO=KK-l 

DO 12 K=1,NV 

B(KK,K»=R(KK,KI/A (KK,KK) 

DO 12 I=l,KHO 
N=KK-I 

DO 11 J=N,KHO 

9 (N,K)=8 «N, K )-A(N, J«-l) *8( J+l, K) 

11 CONTINUE 
BCN,K»=B(N,K)/A(N,N» 

12 CONTINUE 

13 DNATE0=0 

IF (IT.EQ.O) return 
00 3 0 I = 1,JJJ 
DO 3C J=1 , JJJ 


INV 400 
INV 410 
INV 420 
INV 430 
I«V 440 
INV 450 
INV 460 
INV 470 
INV 480 
INV 490 
INV 500 
INV 510 
INV 520 
INV 530 
INV 540 
INV 550 
INV 560 
INV 570 
INV 560 
INV 590 
INV 600 
INV 610 
INV 620 
INV 630 
IMV 640 
INV 650 
INV 660 
INV 670 
INV 680 
IMV 690 
INV 700 
INV 710 
IMV 720 
INV 730 
INV 740 
INV 750 
IMV 760 
INV T70 
IMV 780 
INV 790 


00 



O O O 
O M CM 
OD (BO OD io 


> > > > 
W X 9 W 

M M 
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SUBROUTINE RESPONS RES 

COMMON /NflMES/NN (20J in 9BETANJ11I »DSTftRSll KTPNIll J RES 

1 ,PNF15H 9BETANF<51J »0STARF<51»,PN1F «51» ,BET AN1F!51J* OST R£ S 

2ARlFC51ti sPNC (511 ,BETANC (511 ,OSTARCC 51) ,PN1C(51 > *BETAN1C( 51) *0STAR1RES 
3C(51),APN(11) »ABETANlll) ,A0STAR(11).3PNU1)»B6ETAN(11) ,BDSTAR « 11) » RES 
4PHIN«11) RES 

COMMON /PARAM/MU, EIGEN* C*NCASE*MCASE* N7, NPLOT* RES 

1NHIST*PNFINAL »9NF INAL*DSFINAL»TIME (51) *P(4»4) ,Q C4) RES 

COMMON /LIMITS/PNUC5 1) »PNL (5i)*BETANU(5i)» BETANLC51) *DSTARUC 51) * RES 

1DSTARL(51),PN1U(51). PNl L(51 ) * BETAM UC5 1) .BETANIL (51) *0STAR1U (511 , RES 
2GSTAR1L»51) RES 

COMMON /FINAL/ PSS. BETA SS* 0STARSS*0P,DR*D8ETA * DPMI , A (4, 4) » B( 4) *00£RES 
ILTA RES 

COMPLEX MU(4) »EIGEN(4) *C (5,4) RiS 

OIHENSION TEMP (4, 4), XTEMP(4,4) ,VTeMP(4,4)* R(4*4) RES 

DIMENSION XRK 51) ,XR2(51) ,KR3(51) *XR4(51) R^S 

00 1 1=1,51 R^s 

1 TIMEd ) = FL0AT(I-l)/10o RES 

IF( NN( 9) oEQ,!) CALL LISTER(9) Ri S 

DO 10 1=1,51 RtS 

PNF II ) =C (1,1 )*CEXP(EIGEN m *TIME( I) ) +C(2,1I *CE XP (EIGEN (2 ) *T IME 1 1) ) RES 
1+CC3,1>*CEXP(EIGENC3) »TIMEII) ) +C 1 4 , 1 )*CEXP (EIGEN (4)*TIME (I ) )+-C (5, IRES 


2 ) 

8ETANF (I 
1 (I) ) +C (3 
2C (5,2) 
OSTARFCI 
1 ( I) ) +C(3 
2C(5,4) 
PN1F(I) = 
1*CEXP(EI 
2EIGENI4) 
BETANIFI 
1*CEXP (El 
2EIGEN(4) 
OSTARIFI 
ICEXP (EIG 
2EIGENI4) 
10 CONTINUE 


RES 

)=C( 1,2) *CEXP( EIGEN( 1) ^fTIME(I) H-C (2,2 )*CEXP (E IGEN (2 )* TIMERE S 
,2) »CEXP(EIGFN (3)*TIM£ (I ) )+C (4,2)*CEXP (EIGEN <4)*TIME( I) )♦ RES 

RES 

)=C(l,4)*CEXP(EIGEN(11*TIFE(I) ) *C (2,4>*CEXP (E IGEN (2 )*TIMERE S 
,4)*CEXP(EIGEN(3)*TIME (I) )+G(4,4)*CEXP(EIGEN(4)»TIME(I))+ RES 

RES 

EIGEN(1)*C(1 ,1) ‘»CEXP(EIGEN(l)*TIMEtI) ) +EIGEN ( 2)*C(2, 1 ) RES 
GEN(2)*‘TIME(I) ) + EIGEN(3)*C(3,1 )*CEXP(E IGEN(3 )*TIHE(I) )♦ RES 
*C(4,l)*CEXP(EIGEN(4)*nME(in Ri S 

I)=EIGEN(1)*C( 1 ,2)*CEXP (EIGEN (1)* TIME (1) )+EIGEN (2 )* C ( 2, 2) RES 
GEN (2)*TIME( 1) 1 4.EIGEN( 3 )* C ( 3, 2 ) *CE XP(E IGEN(3 ) *TIME( I) ) ♦ RES 
*C(4,2) *C EXP (EIGEN ( 4) »TIME(I) ) R^S 

I)=EIGEN(1) *C(1 ,4)*CEXP (E I6EN (1 )*TIME ( I ) ) ♦£! GEN ( 2)*C (2, 4) »RES 
EN(2)*TIMEm )-«-EIGEN(3) »C(3,4) *CEXP(E IG EN(3 ) *T IME ( I ) ) * RES 
*C(4,4) *CEXP (EIGEN(4)*TIME(I) ) RES 

RES 


10 

20 

30 

40 

50 

60 

70 

80 

90 

iOO 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

299 

300 
310 
320 
330 
340 
350 
360 
370 
3 80 
390 


00 

40 



(o 

o 


C 

C 

c 



IPtNNtlOl.ea.ll CALL L"STERflOl 

PAYMTERS RECIPE 

AA=Ae3tA(in 
DO 20 J=2,i6 

IF(ABSCA ( jn «GToABS( AA )1 AA^AeSfA'JM 
20 CONTINUE 
AAT=0.1*AA 
FACT=lo 
DO 30 K=l,30 

fact=fact*float<k» 

X = (4.*AAT) ♦♦K^EXP(i*o *AAT>/FACT 
IFCXoLT. 0.001) GO TO 40 
30 NZ=K + 1 
40 CONTINUE 

IF(NN(ll).EQ.l) CALL LlSTER(ll) 

TEMP(1)=TEMP(6)=TEMP(11) =TEMP(16)=1« 

TEMP (2 )=TEMP(3)=TEMP (4) = TEMP ( 5 ) = TEMP ( 7 ) = TEMP ( 8 )=TEMP C 9 )=TEMP (10 ) = | 
TEMP(12> = TEMPC13)=T£MP{ 14) =TEMP(1?) =0. 
P(l,l)=P(2»2)=P(3»3)=P(4,4)=lo 

P (1,2)=P (1,3) = P (1,4)=P (2^ 1)=P(2* 3) = PC2,4) = P( 3, 1) =PI 3 , 2) =P(3,4) = 0. 
P(4,1)=P(4,2)=P(4,3) =0. 

DO 50 I=1,NZ 
DO 51 JJ=1,16 

51 XTEHP( JJ)=0« 0 
DO 55 J=l,4 
DO 55 K=l«4 
DO 55 L=l94 

55 XTEMP(J,K)^XTEMP(J,K )+TEMP(J,L)*A CL v K ) »0. 1 /FLO AT ( I ) 

DO 52 JJ=1,16 

52 TEMP(JJ) =XTEHP( JJ) 

DO 60 11=1,16 

60 P(II)=P( II)«-TEMP(II) 

50 CONTINUE 

IF(NN(12> .EO. 1) CALL LISTERC12) 

TEMPCl ) = TEMP (6)=TEMP ( ill =TEMP( 16) =1. 

TEMP(2) = TEMP(3 ) = TEMPC4) = TEMP (5 ) = TEMP (7 )=TEM»' (6 )=TEMP (9) = TEMP (10 = 1 
TEMP (12) =TEMP(13)=TEMP( lA) = teMP(15) =0* 


RES 

400 

RES 

410 

RES 

420 

RES 

430 

RES 

440 

RES 

450 

RES 

460 

RES 

470 

RES 

4 80 

RES 

490 

RES 

500 

RiS 

510 

RES 

520 

RES 

530 

RES 

540 

RES 

550 

R£S 

560 

RES 

570 

R£S 

3 80 

R£S 

590 

RES 

600 

RES 

610 

RES 

620 

RES 

630 

RES 

640 

RES 

650 

RES 

660 

RES 

670 

RES 

680 

RES 

690 

RES 

70 0 

RES 

710 

RES 

720 

RES 

730 

RES 

740 

RES 

750 

RES 

760 

RES 

770 

RES 

780 



RC1,1} = R« 2,2» = S<3,3J =R«<»,4S = 1. 


RES 

790 


R Ci« 2 »=R fi.3)=R{l94BsRC2«lilsR(2,3)sRS2,4jsR| 3« IB >RC 3 «2 »:^RS3*4B-8« 

RES 

600 


Rt4«H=R<4,2l=Rf4»3>“0. 


RES 

810 


DO 61 I-1»NZ 


RES 

820 


00 62 JJ = l 9 ie 


RES 

830 

62 

XTEMPC JJB-0,0 


RES 

640 


00 65 J=l,4 


RES 

850 


DO 65 K=l,4 


RES 

860 


DO 65 L=l,4 


R£S 

870 

65 

X TEMPI J,KJ=XTEMP fJ,K HTEMPCJ^LJ^A tlL.KI*0, 1/FLOATI I+l > 


RES 

SSO 


DO 67 JJ=1,16 


RES 

890 

er 

TEHPlJJB =XTEMP IJJB 


RES 

900 


DO 69 II-I 9 I 6 


RES 

910 

69 

R(ri>=R(II)+TEMPIIIB 


RES 

920 

61 

CONTINUE 


RE S 

930 


Q(1)=R fl»l )»8 (11 + R 11 ,2 )»B(2 ><^R(l»3)*nC3B+R (1,41*914) 


RES 

940 


Q (?)=R (2,1)*8<1)+R(2,2)*8C2) +Rf 2, 3) *BI3t ♦R(2,4)*8I4) 


RES 

950 


0(3)=R(3 ,1)*8(1)+R(3,2 )*RI2)+R I3,3)*8(3) +R (3,4)*9I4) 


RES 

960 


(H 4 )=R (4,1)*0C1) + RI4,2 )*9(2 ) + R 14 ,3 ) * 8 (3) +R (4,4 ) *9 (4) 


RES 

970 


Q(l)=Qm*0,l 


RES 

980 


Q(2)=QC2)^0.1 


RES 

990 


0 (3)=g (3)*Q.l 


RES 

1000 


Q(4)=Q(4)»0,1 


RES 

1010 


IF(NN(13) oEQ.l) CALL LISTER(13) 


RES 

1020 


XRlIl) =XR3(l)=XR4m=0» 


RES 

1030 


XR2I1)=-D0ELTA/0R 


RES 

1040 


00 70 1=2,51 


RES 

1050 


XRl (I)=P ( 1,1 I*XR1 (I-l) +P(1,2)*XR2 (I-l) +P(1,3) *XR3(I- 1) +P(1,4) *XR4I RES 

10 60 


ll-l) +Q( 1 > 


RES 

1070 


XR2(I)=P (2,1)*XR1(I-1>+P(2,2)*XR2(I-1»+P(2,3)*XR3II-1)+P(2, 

4)*XR4(RES 

1080 


1 I-l) *0(2 ) 


RES 

1090 


XR3(I)=P (3,1) *XR1(I-1)+P(3,2)*XR2 ( I-l )+P (3 ,3 )* XR3 (I-l ) +P (3, 

4)*XR4(RES 

1100 


ll-l ) *0 (3 ) 


RES 

1110 

70 

XR4(I)=P(4,1)*XR1(I-1) +PI4,2)*XR2(I-1) +P ( 4 , 3) *XR3 (I- 1) *P (4, 4)*XR4( RES 

1120 


ll-l) +0(4) 


RES 

1130 


DO 71 1=1,51 


RES 

1140 


PNC(I)=PSS*XR1 (I) 


RES 

1150 


3ETANC (I)=BETASS*XR3(I) 


RES 

1160 

71 

0STARC(I)=DP*XR1(I) +DR*XR2(I) +DBETA*XR 3 ( I) ♦0PHI*XR4 ( I) +DDELTA 

RES 

1170 


iO 

N> 


IF{NN«HfUEQ*l> CALL LISTERJ141 

XRiT=XRifB*A?l*l»*XR2m*A«l,2S*XR3ai)*A«l*3}*XR4a»*An,M*0ll» 
XR2T=XR1 «I I*A «2*1J+XR2 ll »*A (2*2 HXR3 ( 1 1* A «2 13 l♦XR4(Il»•A^2*4l♦B(21l 
XR3T=XRl (IJ*A (3*t 1I + XR2 (I J*A 13* 21 ♦XRSJI I*A( 3*35 *XR4( 1 5 *A( 3*45 ♦B(3I 
XR4T=XR1 C U 1 4 *15 + XR2 ( II *A (4 *2 I +XR3 (I S*A (4*3 5 +XR4(I J* A ( 4*45 *B (45 
PNIC (I l=PSS*XRiT 
BETANICI I5=eETASS*XR3T 

30 DSTAR1C(I»=OP*XR1T4-DR*XR2T+O0ETA*XR3T-i-DPHI*XR4T 

IF(NN(15I eEQoll CALL LISTER(15I 

RETURN 

END 


RES 

RES 

RES 

RES 

RES 

RES 

RES 

RES 

RES 

RES 

RES 

RES 


1180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 



SUBROUTINE MODEL *C,E IGEN.PSS, BETA SS, OSS , DP, OR, DBETA,OPHl . A. 8 , 


iODELTA) 

COMPLEX C?5,4» *EIGENJ<»> 

DIMENSION AJLsAKBIAJ* X tl 6 > , Xl^'Cf 16> ,PC16,16J 
1TIA.45 ,AT14»4J ,HI4» 

real L ,L 1 9 L2 *K * 131 » 1 32 » 1 33 9 134 

DATA X/l8 9 4*0*9 1«9 4*0, 9 1» 9 4*0« 9l*/ 

READ laO 9 V 9 L 9 C 39 GCO 9 PSS 9 BETASS 9 OSS 
print 200,V,L,C3,QC09PSS,BETASS,DSS 

READ 1109lTMAX,EPSl 
PRINT 300,nMAX,EPSl 9 X 
K=C3*QCO 

L1-REALCEIGEN«1H 

L2=REAL«EIGEMC21» 

A3=REAL«EIGENI3H 
B3=AIMAGIEIG£N(3I > 

R11=REAL fC<l9in/PSS 

R12=REAL<C<l92»)/0ETASS 

R13=REAL tC<l93H /PSS 
R14=REAL(Ctl9 4)>/0SS 
R21=RE AL«C{29H?/PSS 
R22=REAL IC(29 211 J /bet ass 
R 23=REALCC<2 93)1 /PSS 
P24=REAL CC(2»4» )/DSS 
R31=REAL (CC39 n »/PSS 
R32=REAL (CC3 92 J > /BET ASS 
R33=REAL (Ct3,3>J/PSS 
R34=RE AHC(3 9 4H /DSS 
I31-AIMA6CC(3,1H/PSS 

I32 = AIMAG(CC3 92H/8ETASS 

I33=AIMAGCC( 3, 3H /PSS 
I34=AIMAGCC (3 94n/0SS 

C1=V*R11 

C2=C1»L/V 

C3=V*R12 

C4=C3*L/V 

C5=V»R13 


PINVtl6»16l*Fll6U 


C6=C5»L/V 
C7= ?R14-K*R12» 


M30 

10 

MOD 

20 

MOO 

30 

MOO 

40 

MOO 

50 

MOO 

60 

M30 

70 

MOO 

80 

MOO 

90 

MOO 

IQO 

HOD 

110 

MOO 

120 

MOO 

130 

HOO 

140 

MOD 

150 

MOO 

160 

MOO 

170 

MOO 

180 

MOO 

190 

MOO 

200 

MOD 

210 

MOO 

220 

MOO 

230 

H30 

240 

MOD 

250 

MOD 

260 

MOO 

270 

MOO 

280 

MOO 

29 0 

MOO 

300 

MOO 

310 

MOO 

320 

MOO 

330 

MOO 

340 

MOO 

350 

MOO 

360 

MOO 

3 70 

H30 

380 

MOO 

390 


C8=»V*F11*L1 

C9sC8*L/V 

01=»C8 

C10 = V*R2i 

C11=C10*L/V 

C12*tf*R22 

C13=C12*L/V 

C14=V*R2 3 

C15=C1 4*L/V 

C16s=CR24-K*R22J 

C17=-V*R21*L2 

C18=C17*L/V 

D2=-C17 

C19=V*R31 

C20=C19*L/V 

C21=V*R32 

C22-C2i*L/V 

C23=V*R33 

C24=C23*L/V 

G25=R34»K*R32 

C26=V*(I3i*83»R31*ft3» 

C27=C26»L/V 

03=-C26 

C28=V*I31 

C29=C28»L/V 

C30 = V* 132 

C31=C30*L/V 

C32=V*I33 

C33=C32*L/V 

C34s^I34-K*I32 

C3S=-V*(R31*33+I31*ft3) 

C36=C35*L/V 

04=-C35 

C37=V*R11 + L*R11*L1 

C38=V»R12+L*R12*L1 

C39=V*R13*L*R13*L1 

C40=V*R11*L1 

C41 = V*R12»L1 

C42=V*R13*L1 


NOD 

400 

H3D 

410 

HDD 

420 

MOD 

430 

NJD 

440 

NOD 

450 

NOD 

460 

NOD 

470 

N3D 

460 

HOD 

490 

HOD 

500 

H30 

510 

NOD 

520 

MOO 

530 

MOD 

540 

HOD 

550 

HOD 

560 

MOO 

570 

MOO 

580 

MOO 

590 

MOD 

600 

MOO 

610 

HOO 

620 

MOO 

630 

MOO 

640 

MOO 

650 

MOO 

660 

MOO 

670 

MOD 

660 

MOO 

690 

MOO 

700 

MOD 

710 

MOO 

720 

MOD 

730 

MOO 

740 

MOD 

750 

MOO 

760 

MOD 

770 

MOO 

780 


D5=R !<>*L l-=-K*R 12*L 1 
C43=V*R21*L*R21*L? 

C 44* V* R2 2*l-'^ R 2 2* L 2 

C45*V*R23*L*R23*«,2 

C46=V*R2i*L2 

C47sV*R22*L2 

C48=V*R2 3*L2 

06=tR24*L2-K*R22*L2 

C49=V*R31+L*R31*A3»L*I31*B3 

C50*V*R32+L*R32*A3»L*I32*B3 

C51=V*R33*L*R33*A3»L*I 33*03 

C52* V*R31* A3 »V* 131*33 

CS3=V*CR32*A3-I32*B3» 

C54=V* R33*A3-T33*B3 > 

Or= R34*A3-I34*B3-K*R32*A3+K*I32*e3 
C55-V* 131 +L*I31*A3<-L*R 31*93 
C56=V* I32+L*R32*03+L*I32*A3 
C57 = V*I33+L*R33*B3+L*I 33* A3 
C58=V*R31*B3 + V*I31*A3 
C59=V* IR 32*93<-I32*A3 > 

C60 = V*<R33*9 3-H33*A3) 

08=R34*B3+I34*A3-K*R32*33-K*I32*A3 

C61=R14-K*R12-V*R12*L1 

C62=-L*R12*L1 

D9=»C62*V/L 

C63=R24»K*R22-V*L2*R22 

C64=-L*L2*R22 

D10 = -C64*V/L 

C65 = R34-K*R32-V*R32*A3+V*I32*B3 
C66=L*<-R32*A3+ 132*33) 

□ 11=-C66*V/L 

C67=I34-K*I32-V*I32*A3-V*R32*B3 

C68*-L*(R32*83+I32*A3) 

O12=-C60*V/L 
CB9=-\/*R13*Ll 
C70=C69*L/V 
013*-C69 
C71=-V*R23*L2 
C72=C71*L/V 
S D14=-C71 


MOO 790 
MOO 800 
810 
MOO 820 
MOD 830 
MOO 840 
HIO 850 
MOO 860 
MOO 870 
H30 880 

MOO 890 
MOO 900 
MOO 910 
MOO 920 
MOO 930 
MOO 940 
M3D 950 
MOO 960 
MOO 970 
M30 980 

MOD 990 
MOD 1000 
H30 1010 
MOO 1020 
MOO 1030 
MOD 1040 
MOO 1050 
MOD 1060 
MOD 1070 
MOO 1080 
MOO 1090 
MOO 1100 
MOD 1110 
MOD 1120 
MOO 1130 
MOO 1140 
H30 1150 
MOO 1160 
MOO 1170 
MOO 1180 


lb 

o» 



C?3=-V*CR33*A3-I33»83 J 

C74=C73*L/tf 

D15*-C73 

C75=-V*«R33^B3^I33*ft3> 

C76=C75*L/V 

D16--C75 

DO h ITER=1»ITMA7 

AA=X MI»X dlOJ»X (2D*X C9» 

AB=Xtl J*X(6>“X(2>*X(5) 

AC=X(3J*XflO) ‘•X(2>*X (111 

aO=X(3J*XC6l-XC2i*X(7l 

AE*X(VI*X(10J=X(2»*X(12» 

aF=XI4»*X(6l"X(2l*X (81 


MOO 

HOO 

MOO 

N3D 

MOO 

HOO 

M30 

HOO 

MOO 

MOO 

HOD 

HOO 

MOO 


HOO 

HOO 

MOO 

HOD 

MOO 

MOO 

MOO 

MOD 

MOO 

MOO 

MOO 


F(l)=Cl*Aa+C2*aBfrC3*AC+C4*aO*C5*aE+C6*AF+Cl*X( lt+C7*X(2H-C3*X(3H- 
1C5*X(4) + C8*X(10H-C9*X(6» -01 

F(2l-C10*Aa4-Cll*AB+G12*AC+C13*AOiC14*AE+C15*AF+C10*X(ll+C16*X«2> + 

lC12*X(3)frC14*X(4)+C17*X (1 0 1 +C18*X (6) -02 
F(3l = C19*AA*C20*AB+C21*aC*C22*AO+C?3*AE*C2 4*AF-»C19*X (1 I +025* X (21 ♦ 
1C21*K(3) +C23*X (41 +C26* X C 1(51 +C27*X (61-03 
F(4»=C28*AA+C29*AB+C30*AC+C31*AO*C32»A€+C33*AF+C28*X C1I+C34*X (21* 
1C30*X(3) +C32*X (4)<-C35*X (10)+C36*X(6)-04 
AA=X(5I*X(10 J-X(6)*X (91 
A9=X(7)*X(10) -X(6»»X(11> 

AC=X(8 l*X (10 l-X (6)*X (12) 

F(5 )=Cl*AA+C3*a3+C5*AC+C37*X(5)+C7»X(6)+C38*X(7) fC39*X (8UC40*K (9) MOO 
1 + C41*X(11H-C42*X(12)-D5 

F (6)=C10*Aa+C12*aB+C14*aC+C4 3»X(5KCl6*X(6HC44*X(7) ♦C45*X(8) ♦C46*M0D 
1 X (9H-C47*X(11 1 ♦C48*X( 121 -05 ^0° 

F f 71 =C19»Aa+C21*A8+C?3*AC+C49»X (5 ) + C25^X (6 »+C50*X(7 I +C5 1*X(8 ) ♦C52*M3D 
1X(9I+C53*X(11I+C54*X(12I-D7 MOO 

F(8I=C28*AA+C30*AB+C 32*AC+C55*X (5) + C34*X (6) +C56*X (71 *057* X(8 M-C58*M00 
IX (9)+C59*X (11» + C60*X (12 I-D8 M3D 

AA = X(9)*X(5)-X(5)*X( 10) MOO 

AB=X(11) *X(6) -X(7)*X (10) MOO 

ac=x ( 12 )*x (6)-X (8)*X (10 ) M30 

F (9)=C2*aA+C4*AB+C6*AC+Cl*X(9H-C61*X(10)+C3*X (IIM-C^^X (121+C62* MOO 


1X16) -09 


MOO 


F(10)=C11*AA+C13*A3«-C15*AC+C1Q»X(9)+C63»X( 10 ) + G12*X ( 11 ) +C14*X (121 +i3Q 


1C64*X(6) -DIO 


MOD 


1190 
1200 
1210 
1220 
1230 
1240 
1250 
1260 
1270 
1280 
1290 
1300 
1310 
1320 
1330 
1340 
1350 
1360 
1370 
138 0 
1390 
1400 
1410 
1420 
1430 
1440 
1450 
1460 
1470 
1480 
1490 
1500 
1510 
1520 
1530 
15 40 
1550 
1560 
1570 





F«llIl-C2 0*flA»C22*flB+C2i»* ftC»C19»X J9 J ♦C65*KI iOB ♦C2l*X« 11 > ♦C23*Xt 12> ♦MOO 
lC6e»X«6S “Oil MOD 

FU2lsC29*flA+C31*AB^C33*ftCfC28*X«9>+C67*X«i9UC30*X«ll>^C32*Xfl2J*MOD 
1C68*XC6>“D12 MOO 

AA=X<13}*XU0J“X«l‘»l*XI9) MOO 

flB=XC13l *X(6»“XC14I) ♦XC5» MOD 

AC=X«15J*X«10J“XJ14I)*X«11> MOD 

AO=X (15» *X<6)“X<14)*X«7» MOO 

AE-XC16»*X(10)“X«14J*X«12l M3D 

AF=X<16»*Xt6>“X<14J*X«8» MOO 

FC13|sCl*AA+C2*AB+C3*AC+C4*ADfrC5»AE^C6*AF+Cl»X tl3J+Cr*X C14UC3* MOD 

1X«15»^C5*X(16J + C69*X ilO J + C70*X(6!“D13 MOO 

Ffl4»=C10*AA+Cll»AB+C12*AC^C13^AD+C14*AE+C15*AF^C10*X813l)+C16* MOO 

lX(14>+C12*X(15l+C14*Kll6)+C71»X(10I+C72*X(6t“D14 HOD 

F (15»=C19*AA+C2C*A8+C21*AC+C22*AD+C23*AE+C24*AF+C19*Xt 13> ♦C25* MOD 

1X(14»+C21*X«15»+C23*X( 16) ♦C73*X 1 10)+C74*X «6)“015 MOO 

FU6l*C2 8*AA+C29* AB+C30*AC>C31*AO*C32*AE4C.33*AF+C2a*Xl 13) ♦CSA* MOO 

lX(14)+C30*Xfl5)+C32*X(16)+C75*X(10)+C76*Xt6)“016 MOD 

DO 5 1=1,255 MOD 

5 PCII=0, MOD 

P (1 ,1)=G1*XC10) ♦C2*X(6)+Ci MOO 

Ptl,2) =-Cl*X *9)-C2*X C5 )“C3»X(11 )-C4*X 87)-C5-»X « 12) -C6* X «8) *07 MOO 

P(1,3)=C3*X (10)+C4*X (6)^C3 MOD 

P(l,4) =C5*XH0) ♦C6*XC6) +05 MOO 

PU,5)=-C2*X(2) MOO 

P(l,6) =C2*X(ll^C4*X<3)+C6*Xt4)+C9 MOD 

P(i,7)=-C4*X<2) MOD 

Ptl,8l=-C5*X{2) MOD 

P81,9)=»C1*X(2) MOO 

Ptl,10) = Cl*Xfl) +C3*X 13)^C5*X (4)+C8 MOO 

P (1,11)=-C3*X C2) MOO 

P81,12) = -C5*XC2) MOO 

P(2,1)=C10*X I10>+C11*X t6)+niO MOO 

P (P 9 2)=-G10*X <9)-Cll*X ( 5)-C12*K(li) -C13*XJ 7)-C14*X(12) -C15*XI8)+ MOO 
1C16 MOO 

P (2,3) =C12*X (10) +C13*X (6)^C12 MOO 

P(2,4)=C14*X (10)+C15*X(6)+C14 MOO 

P(?,5)=-C11*X(2) MOO 

P 82,e)=Cll»X {1)+C13*X(3)+C15»X(4)4C1B MOO 

P (2,7)=-C13*X( 2) MOO 


1580 
1590 
1600 
1610 
1620 
1630 
1640 
1650 
1660 
1670 
1680 
1690 
1700 
1710 
1720 
1730 
1740 
1750 
1760 
1770 
1780 
1790 
18 00 
1810 
1820 
1830 

18 40 
1850 
I860 
1870 
1880 
1890 
1900 
1910 
1920 
1930 

19 40 
1950 
I960 
1970 



00 




Pd2,8Js-C15*Xf2II 

P«2,9J=-C10»X«2I 

P < 2, 10^ =C10*X d til ♦C12*X « 3M-Cli**X «4»*C17 
Pe2,llJ = «-Cl2*X (2» 

P!2,12»s«=*C14*Xf2> 

Pd3,l» =C19*XC10?+C20*XI6>*C19 

PJ3,2)=-C19*XC9]l-C20»XJ5»»C2i*X<llJ-C22*X«7l°C 23*X«12»»C24*ff §>♦ 
1C 25 

P«3,3}=C21*XdlOJ +C22*X C6J+C21 
P «3«4J=C23*X f 10 HC24*X«6HC23 
PI3»5J =-C20*X< 2) 

PC3,6JsC20*X <1J+C22*X «3>+C24*X«41l+C27 

P«3,7>=»C22*Xd2J 

PC3,8J=“C24*X«2) 

Pd3,9J=-C19*XC2> 

P(3, 10} = C19*Xdll)+C21*Xl3)-*-C23*XC41l+C26 

P C3,1U=-C21»X (21 

P«3,12J=-C23*XC2» 

P(4, 11=C28*X<10> + C29*X(6»-^C28 

P«4,2J=»C28*XI9»-C29*X (5»-C30*X tll)»C31*X«ril-C32*X(12»=C33*XC 8»* 

1C34 

Pf 4,3»=C3Q^Xt 10t+C31*XI6M-C30 
P(ifr,4J=C32»X (lQt+C33*X I6KC32 
P(4,5J=-C29*X(2» 

P (4* 6» =C29*X (1 »+C31*X (3) + C33*X (4) +C36 
Pd4,7>=-C31*X(2) 

P«4,8»=-C33*XC2) 

P (4.9J=-C28*X{2D 

Pd4, 10» = C28*XdlJ +C30*X(3»+C32*X(4> + C35 
P<4,11> = -C30*X C2> 

P C4,12)=-C32*X t2 » 

P(5 ,5» =C1*X CIO )+C37 

PC5 »6)--Cl*X(9l»C3*XCll)-C5*XC12) *07 

PC59 71 =C3*X(10H>C38 

PCS ,flJ=C5»X(10l*C 39 

PC5,9»=-C1*X(6) *040 

P(5 9lO» = Cl»XC5»+C3*X (71 frC5*Xte) 

p (5* in=-C3*X( 6J + C41 

P(5,12>=-C5»XC6»+C42 

P(6 95I=C10*X(10» +C43 


noo 

NOO 

H90 

HOO 

NOO 

NOD 

NOO 

NOO 

NOO 

N30 

NOO 

NOD 

NOO 

NOO 

NOO 

HOD 

MOO 

MOO 

N3D 

NOO 

HOD 

MOO 

MOO 

MOO 

N30 

MOD 

MOO 

MOD 

MOO 

MOO 

NOD 

MOO 

MOO 

MOO 

MOO 

MOO 

MOO 

MOO 

MOD 

MOO 


1980 
1990 
2000 
2010 
2020 
203 0 
2040 
2050 
20 60 
2070 
2080 
20 90 
2100 
2110 
2120 
2130 
2140 
2150 
2160 
2170 
2160 
2190 
2200 
2210 
2220 
2230 
2240 
2250 
2260 
2270 
2280 
229 0 
2300 
2310 
2320 
2330 
2340 
2350 
2360 
2370 



P J6,6D s»C10*X8 9!l“Cl2*X« 3L1» = C1***X« 12J4-C16 
P«6,7J sCl2*XU0J+C4<* 

P 86*8J-C14»X 110 HC^5 
P«6*9?-»C10*X«6» ♦C46 

P 86,10 |sC10*X 85»+C12«Xt?>*C14*W 88 J 

P86,11J--C12*XS6KC47 
P86,12l = -C14*X861l *048 
P 8 7 » 5 > =C 19*X 8 1 0 » *049 

P«7,6V=»C19*X89»-C21*XC11J “C23»X812f+C25 

PJ7, 75 =C?1*X(10J +050 
P87,8>=C23*X 810?*C51 
P87,9D=»G19»XC61*C52 
P 8 7, 1011 = C19*X8 5> + C21*X8 7»*C23*X88!) 

P87,11» = -C21*X86»+C53 
P87,121=-C23*X86t*C54 
P88,5)=C28*X810>+C55 

P <8 ,6»=-C28*X 89)-C30*X 8 11» »C32»X812) *C34 

Pf8,7l=C30*X(10»+C56 

P{8,8J=C32»X810)+C57 

P 88,91=“ C28*X(6)-»-C58 

P 8 8, 10)=C28*X8 5)+C30*X 8 7J +C32*X(8» 

P88,11»=-C30»X8 6>+C59 

P88,12I=-C32*XC6H-C60 

P 89, 6) =C2*X89» ♦C4*X( 11) +C6*X8 12) +C62 

P(9,5) = “C2*X 810) 

P (9,7)=-C4»X 810) 

P89,8)=-C6»X810) 

P89,9)=C2*X86M-C1 

P 89,10) = -C2*X 85)-C4*X8 7)-C6*X8 8)+C61 
P89,11)=C4*X86H-C3 
P 89,12) = C6*X (e)*C5 
P810,5)=-C11*X810) 

P810,6)=C11*X89)+C13*X 811)+C15*X C12)+C64 

P 810,8 )=-C15*X 810 

P810,7)=-C13*X810) 

P810,9)=C11*X86) +C10 

P810,1Q)=-C11*X 89)-C13*X87)“C19*X8 8)fC63 
P810,ll) =C13*X86)+C12 
P810,12» =C15*X 86)+C14 


NOO 2380 
HOD 2390 
^3D 2400 
HOD 2410 
MOD 2420 
H30 2430 
MOO 2440 
HOD 2450 
MOO 2460 
hod 2470 
MOO 2480 
MOO 2490 
NOO 2500 
NOO 2510 
NOO 2920 
MOD 2530 
MOD 2540 
NOO 2550 
MOO 2560 
MOO 2570 
MOO 2580 
MOO 2590 
MOO 2600 
MOO 2610 
MOO 2620 
MOO 2630 
MOO 2640 
MOD 2650 
MOO 2660 
MOO 2670 
MOO 2680 
MOO 2690 
MOO 2700 
MOD 2710 
MOD 2720 
MOD 2730 
MOO 2740 
MOO 2 750 
MOO 2760 


o 

o 


PCtl«5»=»C20*Xll0» 

P«ll, 6 I«C 2 D*Xl 9 >*CZ 2 *X!ll>'*‘C 2 «i^K« 12 »»C 66 

P«il»71)--“C22*X«105 

PI11,8}=»C2<»*XC10> 

P«11,9>=C20*X!6I+C19 

P{il,101»--C20*Xf51)-C22«Xf7S»C2i»*XS8>»C65 
Ptll.ll>=C22*X«6)*C21 
P (11,129 =C24*X(69 *023 
P(12,59=-C29*XC109 

PC12,69=C29*X(9»*C31*XI119»C33*X(129»C68 


P(12,79=-C31*X(10I 

Pfl2,89=-C33*X(10) 

P(12,9)=C29*X(6»*C28 

Pfl2,l09=-C29*XC59-C31*X(7»-C33*X(89^C67 
P(12,119=C31*X(69+C30 
P(12,12) =C33*XC6» *C32 
PU3,5J=»C2*X (149 

P (13,6 9 =C2*X( 139 ♦C4*XC 159 ♦C6*X{ 1 69 + C70 
P(13,79=-C4*X(149 
P C13,89=«»C6*X (149 
P(13,99 = -C1*X( 149 

P(13,10 9=C1*X (139 tCl^X (15 9fC5»>(169+C69 

P(13,119=»C3*X(149 

P(13,129=-C5*X(149 


P(13llM=-Cl*x\V9-C2*X(59-C3*X(119-C4*X(79-C5*X(129-C6*X(8 9+C7 

P(13,159=C3*X(109 +C4*X (69+C3 
P (13,16I=C5*X (10 9+C6*X (69+C5 
P (14,59 = -C11*X( 149 

P(14,69 = C11*X (139+C13*X (159*C15*X C16 I + C72 

P(14,79=-C13*X(149 

P(14,99=-CiO*X(149 

P(14,89 = -C15»X (149 

P (14,109 =C10*X(139*C12*X(15J+C14*X(169 *071 

P(14,119=-C12*X(149 

P(14,129=-C14*X(149 

P(14,*149=-C10*X(9 9-C11*X (5^ (119 -C13*>X (7 9 -C14* X (129- C15» X (89 


1C16 


MOO 2770 
MOO 2780 
hod 2790 
HID 2800 
MOO 2810 
HOD 2820 
MOO 2830 
MOO 2840 
HOD 2850 
MOO 2860 
MOO 2870 
MOD 2880 
HOO 2890 
MOO 2900 
MOO 2910 
HOO 2920 
M30 2930 
MOD 2940 
MOO 2 950 
M30 2960 
MOO 2970 
HOO 2988 
HOO 2990 
HOO 3000 
HOD 3010 
MOO 3020 
MOO 3030 
MOO 3040 
MOO 3050 
MOO 3060 
HOO 3070 
HOO 3080 
MOO 30 90 
MOO 3100 
HOO 3110 
MOO 3120 
MOO 3130 
*^MOQ 3140 
MOO 3150 


10 


3 

11 


P«14*1S»=C12*”K«105*C13*X«6J ♦C12 

P«14il6»=C14*X«10J*C15*X«61*C14 

P«15.5 5s»C21II»Xll4J 

Pdl57 6»=C20*X?13? ♦C22*X !15J*C24*X 5168 *€74 

P 515 ,7S=-C22*X 5141 

P515*8B=»C24»X514» 

Pfl5*9J=»C19*X514J 

P515,10 8 = C19*X5138*C21»X515»+C23*X516M€73 

P515,ll) =-C21*X514l 
P515,128-»C23*X5148 

P515*l3»sC19*X510»*C20*XC6J +C19 

P515,148 =-C19*X598»C20*Xt58=C21*X C11J-C22*XC71 
tC25 ^ 

P515,158=C21*X5101+C22*X56I+C21 
P51S*16»=C23*X510>»C24'^X 561 +C23 
P 516, 5»=-“C29*X 5141 

PC16,6) = C29®X513H-C31*X(151 +C33*X516M-C76 
P516 ,71=-C31*X(14) 

P516.8J = “C33*X(14> 

P516,9ls-C28 + X514J 

P516,10»-C28*X513)»C30*X 515 H-C32* X 5161+C75 
P516,11»=»C30»X514» 

P(16»12) =-C32*X514) 

P516 ,13»=C2 8*X C10I+C29*X 561 +C28 

P516,141=-C28*X(9)-G29*X551 -C30*X5111-C31*X57J 

1C34 

P516 ,151=C30*X510H-C31*X561+C30 
P( 16,161 =C32*X5101+C33*X5 61 *032 
2 ITCON=l 

CALL INVR5P,PINV,16» 0» 161 
DO 10 1=1,16 
XlNC5I1=0o 
00 10 J=l,16 

XINC5I1 = XINC 5I1»PIW 51, J1 *F 5J1 
DO 3 1=1,16 

IF 5ADS 5XINC51 1 1 •GToEPSlI ITCON=0 
X5I1=X (I1+XINC5I1 
CONTINUE 

IF( ITCCNsNE. 01 GO TO 1 


HOO 

HOO 

HOD 

HOO 

H90 

HOO 

MOO 

HOO 

M580 

HOO 

»C23*X5121-C24*K58»^H0 0 

HOO 

HOO 

HOO 

HOO 

HOD 

HOO 

HID 

HOO 

HOO 

HJO 

HOO 

HOO 

-C32*X 5121-C33*X581+H30 


HOD 

HOO 

MOO 

HOD 

HOO 

HOD 

HOO 

HOO 

HOO 

M3D 

HOD 

HOD 

H30 

HDD 


3160 
3170 
3180 
3190 
3200 
3210 
3220 
3230 
3240 
3250 
3260 
3270 
3280 
3290 
3300 
3310 
3320 
3330 
3340 
335 0 
3360 
3370 
3380 
3390 
3400 
3410 
3420 
3430 
3440 
3450 
3460 
3470 
3480 
3490 
3500 
3510 
3520 
3530 
3540 



ORIGINAL PAGE IS 
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4 CONTINUE 

PRINT 500*X,XINC 
GO TO 7 

1 PRINT 400*ITERfX 
7 CONTINUE 
DO 6 

00 6 J=l*4 
6 All* J»-X J4*I*»4«-J« 

0P=IV*AC3,1>+L*AC2,1 H 
nR=«V»A«3,2>*L*AC2,2»+V« 

0BETA=fV*AC3, 3M-L*A I 2,3^ +KJ 
0PHI = CV*A«3*4J + L»A(2,4 » S 
00 8 I=lv4 

W<II=lAtI,l>-A «I,2>*0P/0R1*REAL tC <5,1)»/PSS + 

l(A(I,3»-Aa,2i*DBETA/0R»*REAL«C{5,2H/DETASS* 
2 (Afl,4»-fl C 1 , 2 J *DPHI/DR>*REAL(CC5«3n/'PSS* 
3A!l9 2» ♦REALtC <5,4 H /OR /OSS 
6 W(I)=»«(I) 

71=1.”L*A<1,2)/DR 
Z2-»V*A(1,2J /OR 
2 3=-L*A C3,Z>/DR 
Z4-1.»V*A <3,2» /OR 

n (1 ) = <Z 4 *W{ 1 ) - Z2*w<3 n/ (Zi*Z4-Z2*Z3> 
S(3)=<-Z3»WC1)+Z1*W<3) t /(Z1*Z4-Z2*Z3) 

OOEL TA=?y*0(3) ) 

B C2 )-W (2 n-ODELTA^AIZ^ZJ/OR 

B<4» =W(4)+DDELTA*A(4,2 »/0R 

OP=OP*CSS 

OR=DR*OSS 

DBETA=DBETA*OSS 

0PHI=0PHI*DSS 

PRINT 9QO,PSS,BETASS,PSS,OP,DR,08ET A,0PHI 

DDELTA=DDELTA*DSS 

PRINT 8QO,OOELTA 

PRINT 600,(tA(I?J),J=l»4),B(II,I-l,4) 

return 

100 FORMAT <«^FID. 0» 

110 FORMATCI4,F20.0) 


NOD 3550 
HOD 3560 
H30 3570 
HOO 3580 
HOD 35 90 
HOO 3600 
NOD 3610 
MOD 3620 
HOO 3630 
HOO 3640 
HOD 3650 
HOD 3660 
HOO 3670 
HOO 3660 
HOO 3690 
HOO 3700 
MOO 3710 
HOO 3720 
HOO 3730 
HOO 3740 
MOO 3750 
MOO 3760 
MOO 3770 
HOO 3780 
HOO 3790 
HOO 3300 
MOO 3810 
HOO 3820 
HOO 3830 
HOO 3840 
M30 3850 
MOO 3960 
HOO 3870 
HOO 3860 
MOO 3890 
HOO 3900 
H30 3910 
MOO 3920 


200 FORMAT «///,10X,M»HAIRPLANE MODEL WITH SPECIFIED TIME HISTORIES*///MOD 
1 , 10 X 9 29HFLIGMT AND VEHICLE PARAMETERS, //, lOX ,OHAIRSPEEO,F10.1«7M FNOO 
2T/SEC,//,10X,<*1HCG TO PILOT STATION LONGITUDINAL DISTA NCE, F10.2, 3HM0D 
3 FT,//9lOX,39HOIMENSIONAL CONSTANT FOR OSTAR EOUATIOM, F10,<*, 30H CUMOO 
4BIC»FEET/LB=SEC0N0S>=S(JUARE0,//,10X,16H0VNAMIC PRESSURE, FIO, I, 14H LMOD 
5B/FT°SQUARE0,//, lOX , 29 HROLLRATE NORMALIZATION FACTOR,F10.3, ^/, lOX, MOO 
629HSI0ESLiP NORMALIZATION FAC TOR, FIO, 3,//, lOX, 26HDST AR NORMALIZATIMOD 
70N FACTOR, F10,3,/ll NOD 

300 FORMAT «10X,43HMAXIMUM NUMBER OF NEMTON-RAPHSON ITERA TI ONS, 18, / , 10XH3O 
1 ,17HC0NVERGENCE LIMIT, F20.5MOO 

2, X ,16F8.4,//J NOD 

400 FORMAT *///,10X,25HITERATIONS TO CCNVERGENCE,I0 ,//,10k,12HN“R SOLUTMOO 
1I0N,//,4X,16F8.4,///J NDO 

EOO FORMAT flOX,29HN-R SOLUTION DID NOT CON VERGE,/, lOX, 19MLAST TRIAL SOMOO 
1LUTION,/,4X,16F8,4,/,10K,24HLAST SOLUTION A 0 JUSTMENT ,/,4X,16F6 .4? M3D 
600 FORMATtlOX,lEHTHE A MATRIX IS ,59X ,15HTHE B VECTOR IS // ,4 (lOX ,4E16. MOD 


16 , 10 X,E 16 . 6 //n 

700 FORMAT(//,10X,20HMO0IFIED E IGENMA TRIX, // ,4 ClOX ,4E20 o 6 , //U MOD 

800 FORMAT? lOX , 14H(0ST AR J OELT AA= ,F 20 * 6 , ///? MOD 

900 FORMAT l//,10X,19HDISTRI8UTION MATRIX,//, 10X,F20,4,//,5 OX,F20, 4, MOO 
1//,70X,F20. 4,//,10X, 4E2 0.6,//» N30 

1 300 FORMAT?10X,4E?0,6» NOO 

END nod 


3930 
3940 
3950 
3960 
3970 
3980 
3 990 
4000 
4010 
4020 
4030 
4040 
4050 
40 60 
4070 
4080 
40 90 
4100 
4110 
4120 
4130 
4140 
4150 



APPENDIX C 

Program AANDB Sample Output 



! 


1 


j 


I 


DATES 12/01/75 


{LIST 

PARAMETER 

table i 

1 1 

111 

1111 

1 1 

CASE NUHBER 1 PLOT CODE 

INPUT DATA 

0 

TIME i 

RESPONSE CODE 

1 

PN» 

0.000 

• 620 

1.020 

1.0 20 

1.005 

1.0 30 

TINE, 

0.00 

.50 

1,00 

1.50 

2 . or 

2.50 

BETAN» 

0.000 

.160 

• 530 

.9 30 

,980 

.800 

TIME, 

0.00 

.50 

1.0 0 

1. 50 

2.00 

2,50 

OSTAR, 

0. 000 

.376 

.676 

.951 

1. 26 0 

1.610 

TIME, 

0,00 

. .50 

1.00 

1.50 

2.0C 

2.50 


1 1 

111 

1 1 



1,070 

3.00 

1.090 

3.50 

1.075 

4.00 

1.045 

4.50 

1.025 

5,00 

• 650 
3.00 

• 690 
3.50 

. 890 
4.00 

1.090 

4.50 

1.150 
5. 00 

1.980 

3.00 

2.3A0 

3.50 

2.650 
4. 00 

2.950 

4.50 

3.270 

5.00 


105 



REAL 


IHA6 


REAL 


I MAG 


HUI» 

£!• 


REAL 

IHAG 

REAL 

IMAG 

• 31766 2E«’00 

• 749635E+00 

• 317662E + 00 

-,7 49635E*00 

-.4112 

2,3400 

-•41J.2 

-2.3400 


.159<»3€E + 01 Q. 

•9330 Q»OOQO 


.279948Ef00 Oo 
=2.5470 0.0000 


REAL IMAG 

MU, .932945E + 00 0. 

El, =.i388 0,0000 


REAL IMAC- 

=a515220£+0l Oo 

3.2783 6.2832 


REAL 

•497955E+00 

-.1528 


IMAG 

0 78125 6E*00 
2.00&7 


REAL 

• 497955E + 00 
=• 1528 


IHAG 

-• 781256E+00 
-2,0067 


original Page is 

OP POOR QUAUi^i 


106 



REAL 

IMAG 

PEAL 

IHAG 

HU« 

• 103636E«Qi 

0. 

,3926946+00 

0. 

El. 

• 071 <» 

0.0000 

“1.8694 

0.0000 


REAL 

IKAG 

REAL 

IMAG 


• ?56iiS£4^00 

•105674E^01 

.2561186 + 00 

-.105674E+01 


• 1675 

S , 6660 

.1675 

•2.6660 


107 


MU VALUES SPECIFIED 


REAL 

IN AG 

REAL 

IHAG 

MU« •3005i7E*OD 

0. 

,998453E*-00 

0. 

EIGEN, 

0*0000 

». 0031 

0.0 000 

REAL 

IMAG 

REAL 

IMAG 

• 4S1S62E«0Q 

».756023E+00 

• 4515e2E«-QQ 

•756023E+00 

-,2$43 

-2.0&48 

“•2543 

2*0648 


THE GENERATED PHI DATA IS, 


PHIN, 

TINE, 

0,000 

0*00 

*247 
• 50 

*721 
1* 0 0 

1*235 
1* 50 

1.740 

2.00 

2*247 

2.50 



2*772 

3,00 

3.313 

3,50 

3.856 

4.00 

4*386 
4* 5 0 

4*905 

5*00 


^ ^^OOR qualj ^ 


j08 


TIME RESPONSE COEFFICIENTS 



REAL 

INAG 

REAL 

INAG 

REAL • 

1NA6 

PNv 


» • 000 0 

2,1496 

.0000 

• 015S 

• 0414 

BET AN , 


°«aooo 

«32,9159 

,0000 

-.1968 

-•1282 

PH IN, 

• 443S 

=.nooo 

»345,4494 

• 0000 

-.0190 

• 0034 

tJSTAR, 

».eC86 

OOQC 

-215,988r 

,0000 

.0 1T9 

• 0368 


REAL 

IHAG 

REAL 

INAG 

.0155 

-.0414 

V 

—1 . 0 695 

- .OCOO 

-.1983 

• 12 82 

33.4334 

-.0000 

-.0190 

-.0034 

345.0439 

- o 0 too 

.0179 

-.0 368 

215.9615 

— • 0 000 



FITTING RESULTS 



REAL 

IHAG 

REAL 

IHAG 

REAL 

IHAG 

NU« 

• 300S 

0,0800 

• 9985 

0,9000 

•4516 

=,7560 

EIGEN, 

=2,4045 

0.0000 

=•0031 

0.0000 

=•2543 

=2.0648 

COEFFICIENTS 






PN, 

=1,1111 

-•0000 

2,1496 

.0000 

• 0155 

.0414 

BETAN , 

=•1199 

=,0000 

-32,9159 

• 0000 

-•1980 

-.1282 

PHIN, 

• 4435 

-.0000 

-3 45, 4494 

• 0000 

-.0 190 

. 0034 

OSTAR, 

-.0086 

-•0000 

-215.9887 

.0000 

.0179 

.0368 


RE4L IMAG t^AL IHAG 

•^S16 .7560 

^•2543 2,06*8 


.0166 

».1988 

-.0190 


0414 

-1 

12S2 

33 

0034 

345 

03ES 

215 


• 0 696 -.0 000 

•4334 -.8000 

•0 439 -.0880 

•9616 -.0000 


• 0179 



Kimnd Hoo<j 
K a^)vd Tvmomo 


ftIRPUNE MODEL MITH SPECIFIED TIME HISTORIES 

FLIGHT ftNO VEHICLE PARAMETERS 
AIRSPEED 61?. 2 FT/SEC 

CG TO PILOT STATION LONGITUDINAL DISTANCE 22.24 FT 

DIMENSIONAL CONSTANT FOR OSTAR EQUATION -.3190 CUBIC-FEET/L8-SEC0NDS-SQUARE0 
DYNAMIC PRESSURE 331.8 LB/FT-SOUAREO 

ROLLRATE NOR MALI2ATI0N FACTOR .500 

SIDESLIP' TIDRMACI7ATrON FACTOR 19. 000 

DSTAR NORMALIZATION FACTOR .010 

MAXIMUM NUMBER OF NEWTON-RA PHSON ITERATIONS 50 

CONVERGENCE LIMIT .000010 

X l.GOOO 0.0000 0,0000 0.0000 0,0000 1.0000 0.0000 0,0000 0.0000 0.0000 

l.QQQO O.OOQQ 0.0000 0.0000 0.0000 1.0000 


ITERATIONS TO CONVERGENCE 


■-■■■■V ■ ■ . 1 ... 




I . ... .. - 1-.^ - J 



U2 


DISTRIBUTION MATRIX 


& 

o 

o 

o 

e 

♦ 



• 



CM 




%o 


CM 


w 




o 

vi 

e 

o 

o 

4 

o 

LsJ 

m 


o 


•-» 

sn 








0 
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THE A HATRIX IS 


-•235877E+01 
-.53A355 E-Ol 
• 253 874 E- 01 
•959579 E^OO 


•770fl32E*00 

««36869(}Ef00 

••10Q515E*Oi 

».207575E»01 


-•10939 0E+C2 
• 384788E*01 
-•201315E+C0 
• 278157e*C0 


-•2556205-02 
-•20 6077E-03 

• 532788E-C1 

• 262703E-02 


THE 0 VECTOR IS 
•5SlB00E»0t 
-•527I09E-01 
-•21S179E+00 
•495S07E-01 
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RESPONSE ENVELOPES 


TINE 

0. flO 

• iO 

• ED 

• 30 

• <»0 
• 60 
.60 

• 70 

• ao 

.90 

1 . 00 

1 . 10 

1 . EO 
1*30 

1.40 

1 . 50 
1.60 

1.70 
1«80 
1*90 
2.00 

2 . 10 

2 . 20 

2 . 30 

2.40 

2.50 
2 , 60 

2.70 
2.60 

2.90 
3.00 

3 . 10 
3.20 

3 . 30 

3 . 40 

3.50 

3 . 60 

3.70 
3.80 

3.90 

4 . 00 
4.10 

4 . 20 

4 . 30 
4 . 40 

4.50 
4 . 60 


PN BETAN OSTAR 


0.000 

0,000 

2.000 

• 333 

• 010 

2.000 

• 667 

.027 

2.000 

1.000 

« 065 

2.000 

1.043 

.107 

2.000 

1.064 

• 160 

2.000 

1.076 

. 227 

2.000 

1.080 

• 300 

2.000 

1.080 

. 360 

2.000 

1.078 

• 409 

2.000 

1.076 

.450 

2.000 

1.074 

.479 

2.000 

1.0 72 

.507 

2.000 

1.070 

.530 

2.000 

1.069 

.556 

2.000 

1.067 

.579 

2.000 

1.065 

.600 

2.000 

1.063 

.620 

2.000 

1.061 

.639 

2.000 

1.059 

.656 

2.000 

1.057 

.677 

2.000 

1.055 

.692 

2.100 

1.053 

. 70 8 

2.200 

1.051 

.722 

2.300 

1.050 

.733 

2.400 

1.046 

.748 

2.500 

1.046 

.762 

2.600 

1.044 

.778 

2 . 700 

1.042 

.790 

2.800 

1 . 040 

. 80 4 

2.900 

1.038 

.816 

3.000 

1.036 

.829 

3.100 

1.034 

. 8 • . 0 

3.200 

1.0 32 

.852 

3.300 

1.030 

.861 

3.400 

1.029 

.871 

3.500 

1.027 

. 680 

3.600 

1.025 

• 888 

3.700 

1.023 

.896 

3,800 

1.021 

. 90 8 

3.900 

1.019 

.915 

4.000 

1.017 

.920 

4.100 

1.015 

.929 

4.200 

1.013 

.93 8 

4.300 

1.011 

.947 

4.40 0 

1.010 

.956 

4.500 

1.006 

.965 

4.600 


• 2.0 00 

' 4.000 

- 4 . 0 00 

- 2,000 

4.000 

- 4,000 

• 2.0 00 

4.000 

* 4. 000 

• 2,000 

4.000 

• 4.000 

- 2.0 00 

4.000 

- 4 , 000 

• 2,000 

4.000 

- 4.000 

- 2.0 00 

4.000 

- 4.000 

- 2.000 

4 . 000 

- 4 . 000 

• 2.0 00 

4.000 

- 4 . 0 D 0 

- 2 . 0 00 

4,000 

- 4 . 000 

- 2.0 00 

4.000 

- 4.000 

- 2,0 00 

4.000 

- 4.000 

- 2.000 

4.000 

- 4 . 000 

- 2 . 000 

4.000 

- 4 .0 00 

- 2.000 

4.000 

- 4.0 00 

- 2.000 

4.000 

- 4,000 

- 2.000 

4.000 

- 4.000 

- 2.0 00 

4.000 

- 4.000 

- 2.0 00 

4.000 

- 4 . 0 00 

- 2.0 00 

4.000 

- 4 . OOS 

- 2.000 

4.000 

- 4.000 

- 2.100 

4. 200 

- 4.200 

- 2.200 

4.400 

- 4.400 

- 2,3 CO 

4.600 

- 4.600 

- 2.400 

4,800 

- 4,800 

- 2.500 

5.000 

- 5.0 00 

- 2.6 00 

5.200 

- 5.200 

- 2.700 

5.400 

- 5,400 

= 2,800 

5.600 

- 5.600 

- 2.900 

5.800 

- 5.8 00 

- 3.0 00 

6,000 

- 6.000 

- 3.100 

6,200 

- 6,200 

= 3.200 

6,400 

- 6.400 

- 3.300 

6,600 

- 6.600 

= 3.400 

6.800 

- 6.800 

- 3,5 00 

7.000 

- 7.0 00 

- 3.600 

7.200 

- 7.200 

- 3.700 

7.400 

- 7.400 

- 3.800 

7.600 

- 7,6 00 

- 3.9 00 

7,800 

- 7.800 

- 4 . 0 00 

8.000 

- 8.000 

- 4.100 

8.200 

- 8.200 

- 4. 200 

8,400 

- 8 . 4 G 0 

- 4,3 00 

8.600 

- 8 . 6 G 0 

- 4 . 400 

3.800 

- 8,800 

- 4,5 00 

9.000 

- 9.000 

- 4.600 

9 .200 

- 9.200 
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J 


I 


1 . .. • 


4.70 1.006 

4*60 1.004 

4.90 1. 002 

$•00 1.000 


.974 
.983 
.992 
1 . 000 


4.700 

4.600 

4.900 

5.000 


«»4. 7 00 
»4.8C0 
-4.9 00 
-5.0 00 


9o40() 

9.600 

9.800 

10.000 


-9.400 
-9. 600 
- 9.8 00 
- 10. 000 
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1 


J 


I 


PND 07 


4.000 

• 806 

4 . 000 

• 711 

i.€86 

.255 

1.122 

••224 

.992 

-.297 

• 908 

-.323 

. 843 

-.327 

• 797 

-.323 

• 724 

••319 

■ •sei 

-.304 

• 613 

-.289 

• 996 

-.274 

• 9 i 0 

-.251 

• 464 

-.228 

• 421 

-.205 

• 379 

-.183 , 

• 349 

-.148 

• 31 B 

-.129 

.295 

-.125 

.272 

-.118 

.257 

-.125 

• 234 

-.123 

.222 

-.121 

• 211 

-.120 

.199 

-.118 

.192 

-.116 

• 188 

-.114 

.164 

-.112 

.180 

-.110 

.176 

-.109 

.172 

-.107 

• 168 

-.105 

.165 

-.103 

• 169 

-.102 

.165 

-.100 

• 165 

-.098 

.157 

-.096 

.142 

-.095 

. 119 

-.0 93 

• 107 

-.091 

• 096 

-•069 

.080 

-•088 

.073 

-.086 

• 069 

-.084 

• 061 

-.082 

.054 

-.0 81 

• 046 

-.079 

• 046 

-.077 

.046 

-.075 


OETRNDQT 


6.000 

- 6 .000 

6.000 

- 6.000 

4.800 

—4 #8 00 

4.150 

- 4,150 

3,720 

- 3.720 

3.280 

- 3.280 

2.910 

- 2.910 

2.610 

- 2.610 

2.370 

- 2.370 

2 . 100 

- 2.100 

1.950 

- 1.950 

1.780 

- 1.780 

1.610 

- 1.610 

1.490 

- 1,490 

1.350 

= 1.350 

1.260 

- 1,260 

1.170 

- 1,170 

1.120 

- 1,120 

1.070 

- 1.070 

1 . 050 

- 1,050 

1.000 

= 1.000 

1.000 

- 1.000 

1 . 000 

- 1.000 

1.000 

- 1.000 

1.000 

- 1.000 

1.000 

-loOOO 

1.000 

- 1,000 

1,000 

- 1 . 000 

1,000 

- 1.000 

1.000 

- 1. 000 

l.OOD 

- 1,000 

l.OOD 

- 1.000 

1.000 

- 1.000 

1.000 

-1 .000 

1.000 

- 1,000 

l.OOC 

- l.OCO 

1.000 

- 1.000 

1.000 

- 1.000 

1.000 

- 1.000 

1.000 

- 1.0 00 

1.000 

- 1.000 

1.000 

= 1.000 

1.000 

- 1,000 

l.OOD 

- 1.000 

1.000 

-loOOO 

1.000 

- 1.000 

1.000 

- 1.000 

1,000 

- 1. 000 

1.000 

- 1.000 


OSTfiRDOT 


12.000 

• 12.000 

12.000 

- 12.000 

9 . 66 C 

- 9.600 

6.300 

- 8.300 

7.440 

- 7,440 

6.560 

- 6.5 60 

5.820 

- 5.8 20 

5.220 

- 5.220 

4.740 

- 4.740 

4.200 

- 4 . 2 G 0 

3.900 

- 3.9 00 

3.560 

- 3.560 

3.220 

- 3.220 

2.980 

- 2,980 

2.700 

- 2.700 

2.520 

- 2.520 

2.340 

- 2.340 

2.240 

- 2.240 

2.140 

- 2.140 

2.100 

- 2.100 

2,000 

- 2.0 00 

2.000 

- 2.000 

2.000 

- 2.0 00 

2.000 

- 2.000 

2.000 

- 2.000 

2.000 

- 2.0 00 

2.000 

- 2.000 

2.000 

- 2.000 

2,000 

- 2.000 

2.000 

- 2.0 00 

2.000 

- 2,0 00 

2.000 

- 2.000 

2.000 

- 2.0 00 

2.000 

- 2,000 

2.0 00 

- 2,000 

2.000 

- 2.000 

2,000 

- 2,000 

2.000 

- 2.000 

2.000 

- 2.0 00 

2,000 

- 2.000 

2.000 

- 2 . 0 0 0 

2.000 

- 2.000 

2.000 

- 2.0 00 

2,000 

- 2.000 

2.000 

- 2,0 00 

2,000 

- 2.000 

2.000 

- 2,000 

2.000 

- 2.000 

2.000 

- 2.0 00 
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FITTED TIME RESPONSES 


TIME 

PN 

0. 

• 000 

• 10 

• 252 

• 20 

• 450 

• 30 

• 606 

• 40 

.727 

• 50 

• 019 

• 60 

.889 

• 70 

• 941 

• 80 

.978 

• 90 

1.003 

t. 00 

1.018 

1«10 

1.027 

1.20 

1.030 

1*30 

1.029 

• 1.40 

1.026 

1.50 

1*021 

1.60 

1.017 

1.70 

1.012 

1.80 

1.009 

1.90 

1.006 

2.00 

1.006 

2. 10 

1.007 

2*20 

1.010 

2. 30 

1.015 

2,40 

2.021 

2.50 

1.028 

2.60 

1*037 

2*70 

1.045 

2« 80 

1.054 

2,90 

1.063 

3. 00 

1,070 

3. ID 

1,077 

3.20 

1.083 

3.30 

1.087 

3,40 

1.09 0 

3. 50 

1.092 

3.60 

1.091 

3.70 

1.090 

3.80 

1.086 

3.90 

1.082 

4. 00 

1.077 

4. 10 

1.071 

4. 20 

1.064 

4.30 

1.057 

4*40 

1.050 

4.50 

1. 044 

4.6 0 

1.038 


0ETAN OSTAR PNOOT 9ETANDOT 


».000 

-.000 

. 003 

. 082 

• 020 

• 160 

.052 

• 236 

• 099 

.308 

. 159 

.376 

.231 

.441 

.312 

.504 

. 399 

. 563 

.489 

• 621 

,579 

.677 

» 666 

.732 

• 747 

.786 

,819 

.840 

. 881 

.895 

.930 

. 951 

. 9 66 

1.0 09 

. 989 

1.066 

.998 

1.129 

. 995 

1.192 

.980 

1.257 

, 955 

1.325 

.923 

1. 394 

• 865 

1.465 

. 844 

1.538 

.802 

1.611 

.761 

1.686 

.723 

1. 761 

.691 

1 .835 

. 665 

1. 910 

.648 

1.983 

. 639 

2.0 56 

.639 

2.128 

. 648 

2,198 

.666 

2,267 

.691 

2.334 

.724 

2, 400 

.762 

2,464 

0 804 

2.527 

.848 

2.590 

0 894 

2.651 

. 939 

2.712 

, 982 

2,772 

1. 0 22 

2, 832 

1.057 

2.893 

1. 0 88 

2.953 

1.112 

3.015 


2.828 

-.038 

2.232 

.097 

1.754 

• 245 

1.369 

.395 

1.058 

.538 

• 806 

.665 

.601 

, 769 

.436 

.846 

• 3 0 4 

0 893 

.199 

.907 

,117 

,690 

.055 

. 843 

, 009 

.769 

-.021 

. 671 

-.040 

.556 

-.048 

. 428 

-.048 

. 294 

“ « 0 41 

.158 

-.0 29 

. 028 

-.013 

093 

.004 

199 

.022 

-» 288 

• 039 

-.355 

.055 

400 

.068 

A21 

.078 

-. 419 

. 085 

-. 396 

• 068 

-. 352 

, 087 

29 2 

.0 82 

-.218 

. 075 

134 

.0 64 

-.045 

.051 

. 046 

.037 

.135 

.021 

. 218 

.005 

.292 

•.Oil 

,355 

-.025 

, 403 

“.038 

.437 

“.049 

. 454 

* 6 058 

.456 

• # 064 

. 443 

*.067 

0 416 

•.068 

.377 

>,066 

, 329 

> . 0 62 

,273 

-.055 

, 213 


OSTARDOT 

.832 
.803 
.770 
.736 
.70 2 
.668 
• 637 
.609 
.585 
. 566 
.55 3 
. 545 
.542 
,545 
.553 
.566 
. 592 
.601 
.621 
.643 
. 664 
.694 
.703 
. 719 
.731 
,741 
.746 
.748 
.746 
.741 
.733 
.722 
.709 
.695 
.680 
.665 
.651 
.638 
.626 
.617 
.609 
,6 05 
.603 
.603 
. 606 
.611 
.617 
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1 


) 



k, 70 
^•80 
4,90 
9.00 


i«Q33 

1,029 

1«62& 

1,024 


1,130 
1,142 
1,148 
1, 149 


3,077 

3.140 

3,204 

3.269 


»,047 
°,037 
-.0 27 
»,016 


.151 ,62 

, 090 ,63 

,032 ,644 

-,019 ,654 
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PaVWTERS RECIPE MUM8ER IS* 


19 


DIFFERENCE EQUATION P HATRI«* 


• 768 

.117 

-.9i»4 

-•00<» 

• 946 

.374 

• 003 

-.097 

• 960 

• 085 

• 002 

-.021 


DIFFERENCE EQUATION Q VECTOR* ,528 

-.011 

-.020 

.031 


-.0 C 3 
• 0 01 
• 0 09 
1.0 00 


I 


i 


INTEGRATED 

TIME 

0 * 00 
• 10 
• 20 

• 30 

• ()0 
• 60 
• 60 

• 70 
.80 

• 90 
!• 00 
1 . 10 
1.20 
&• 30 

1.40 
1.60 
1.60 
1 . 70 
1.60 

1.90 
2 . DO 
2 . 10 
2 « 20 

2.30 

2 . i «0 
2 « SO 
2.60 

2.70 
2 . 60 

2.90 
3.00 

3 . 10 

3.20 
3*30 

3.40 

3 . 60 

3.60 

3.70 
3*60 

3.90 

4 . 00 
4 . 10 

4.20 

4.30 

4.40 

4.60 
4 * 60 


..;L 


RESPONSE 


PN 

0.000 
.252 
.450 
.606 
.727 
• 819 
, 889 
.941 
.978 
1.003 
1.016 

1.027 
1.030 
1.029 . 

1.0 26 

1.0 21 
1.017 
1.012 

1.0 09 

1.0 06 

1.0 06 
1.007 
1.010 
1.015 

1.0 21 

1.028 

1.037 
1.045 
1.054 

1.0 63 

1.070 

1.0 77 
1.083 
1 . 08 T 

1.090 

1.0 92 

1.091 
1.090 
1.086 
1,082 

1.0 77 

1.071 

1.0 64 
1,057 
1.050 
1.044 

1.038 


BET AN 

Q.OQO 

• 003 

• 020 
.052 
.0 99 
.159 
.231 
,312 
.399 
.489 
.579 
. 666 
.747 
.819 
.881 
.930 
.966 
.989 
.998 
.995 
.980 
.955 
.923 
,885 
,844 
.802 
,761 
.723 
.691 
.665 
.648 
.639 
,639 
,648 
.666 
.691 
.724 
.762 
.804 
.848 
. 894 
,939 
.982 

1,022 

1.057 

1.088 

1.112 


DSTAR 

-.000 
• 082 
.160 
, 236 
,3 08 
,3 76 
.441 
.5 04 
.563 
,521 
.677 
,732 
.786 
.8 40 
0 895 
,951 
1.009 
1 . 068 
1.129 
1.192 
1.257 
1,325 
1,3 94 
1.465 
1,538 
. 1 . 611 
1,686 

1 , 761 
1.335 
1,910 
1,983 
2,056 

2 . 128 
2.198 
2.267 
2.334 
2,400 
2,464 
2,527 
2.590 
2,651 
2,712 
2.772 
2 , 8 32 
2.893 
2.953 
3.015 


^GOR i 
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4.7Q 

1.033 

1.130 

3.0 77 

4. SO 

1.0 29 

1.142 

3.140 


1.0 26 

1.148 

3.204 

5*00 

1.0 2<» 

1.149 

3. 269 
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1 


I 




.1 J. 


FIRST OER WAT IVES OF 

TIME 

PNDOT 

0.00 

2.328 

• IQ 

2.232 

• 20 

1.754 

• 30 

1.369 

• 40 

1.058 

• 50 

• 8 06 

• 60 

.6 01 

• 70 

• 436 

• 60 

• 3 04 

• 90 

.199 

i. 00 

.117 

!• 10 

.0 55 

!• 20 

.0 09 

!• 30 

“•0 21 

1.40 

-.0 40 

!• 50 

■»* 0 48 

1.60 

-.048 

1.70 

-.041 

i. 30 

-.029 

1.90 

-.013 

2. 00 

• 004 

2« 10 

• 0 22 

2. 20 

.039 

2.30 

.0 55 

2.40 

• 068 

2.50 

.078 

2.60 

.085 

2.70 

.0 88 

2.80 

.087 

2. 90 

• 0 32 

3. 00 

.0 75 

3. 10 

• 0 64 

3.20 

• 0 51 

3. 30 

.037 

3. 40 

.021 

3.50 

.005 

3.60 

-.0 11 

3.70 

-.0 25 

3.60 

-.038 

3.90 

-.049 

4. 00 

-.058 

4. 10 

-.0 64 

4.20 

-.067 

4. 30 

-.068 

4. 40 

-.066 

4.60 

-.0 62 

4.60 

-.055 


INTEGRATED RESPONSES 


BETANO OT 

DSTARDOT 

-.038 

,832 

.097 

.803 

.2 45 

.770 

.395 

.736 

.538 

.702 

.665 

. 666 

.769 

,637 

• 846 

.609 

.893 

.585 

.907 

.566 

,890 

.553 

.843 

.545 

.769 

.542 

,671 

.545 

.556 

.553 

.420 

.566 

.294 

6 3 82 

.158 

.601 

.028 

.621 

-.093 

« 643 

-.199 

• 664 

-.288 

.684 

-.355 

.7 03 

-.400 

.719 

-.421 

,731 

-.419 

,741 

-.396 

.746 

-.352 

, 748 

-.292 

.746 

-.218 

.741 

-.134 

, 733 

-.045 

,722 

.046 

.709 

.135 

.695 

.218 

.6 80 

.292 

.665 

.355 

.651 

.403 

.638 

.437 

. 626 

.454 

.617 

,456 

.609 

.443 

. 605 

.416 

.603 

.3 77 

, 603 

.329 

.606 

.273 

.611 

,213 

.617 


424 



4 t 70 

“.0 47 

ol 51 

• 629 

4«80 

*•037 

.090 

• 634 

4.90 

-•027 

,032 

.644 

5«00 

-•0 16 

= ,019 

0 694 
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APPENDIX D 


Listing of Program AANDB 














LEHIGH UNIVERSITY 
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